You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
317 lines
11 KiB
317 lines
11 KiB
#!/usr/bin/env julia
|
|
|
|
using DataFrames
|
|
using Dates
|
|
using CSV
|
|
using EzXML
|
|
using Images
|
|
using Colors
|
|
#using ImageDraw
|
|
using CairoMakie
|
|
using Printf
|
|
using TimeZones
|
|
using PlotUtils
|
|
#using Profile
|
|
|
|
struct StationInfo
|
|
id::String
|
|
name::String
|
|
num::Int64
|
|
lat::Float32
|
|
lon::Float32
|
|
end
|
|
|
|
function get_station_info(station_id)
|
|
stations = readxml("activestations.xml")
|
|
station_node = findfirst("/stations/station[@id=\"$station_id\"]", stations)
|
|
|
|
#<station id="dukn7" lat="36.184" lon="-75.746" elev="7.7"
|
|
# name=" 8651370 - Duck Pier, NC " owner="NOS" pgm="NOS/CO-OPS" type="fixed"
|
|
# met="y" currents="n" waterquality="n" dart="n"/>
|
|
name = strip(station_node["name"])
|
|
num, name = split(name, " - ")
|
|
lat, lon = parse.(Float32, [station_node["lat"], station_node["lon"]])
|
|
return StationInfo(station_id, name, parse(Int, num), lat, lon)
|
|
end
|
|
|
|
function get_station_lat_lon(station_id)
|
|
stations = readxml("activestations.xml")
|
|
station_node = findfirst("/stations/station[@id=\"$station_id\"]", stations)
|
|
return parse.(Float32, [station_node["lat"], station_node["lon"]])
|
|
end
|
|
|
|
tryparsem(T, str) = something(tryparse(T, str), missing)
|
|
|
|
function get_tile_xy_fraction(lat, lon, zoom)
|
|
n = 2^zoom
|
|
lat_rad = lat * pi / 180
|
|
xtile = (lon + 180) / 360 * n
|
|
ytile = (1 - log(tan(lat_rad) + 1 / cos(lat_rad)) / pi) / 2 * n
|
|
return xtile, ytile
|
|
end
|
|
|
|
function get_tile_xy(lat, lon, zoom)
|
|
xf, yf = get_tile_xy_fraction(lat, lon, zoom)
|
|
return (Int(floor(xf)), Int(floor(yf)))
|
|
end
|
|
|
|
function get_tile(x, y, zoom)
|
|
url = "https://tile.openstreetmap.org/$zoom/$x/$y.png"
|
|
tmp_path = download(url)
|
|
tile = load(tmp_path)
|
|
return tile
|
|
end
|
|
|
|
function get_tile_by_lat_lon(lat, lon, zoom)
|
|
xtile, ytile = get_tile_xy(lat, lon, zoom)
|
|
return get_tile(xtile, ytile, zoom)
|
|
end
|
|
|
|
function get_water_levels(station_number, start_date, end_date;
|
|
timezone=tz"America/New_York", prediction=false)
|
|
start_s = Dates.format(start_date, "yyyymmdd")
|
|
end_s = Dates.format(end_date, "yyyymmdd")
|
|
if prediction
|
|
prod = "predictions"
|
|
else
|
|
prod = "water_level"
|
|
end
|
|
url = "https://api.tidesandcurrents.noaa.gov/api/prod/datagetter?product=$(prod)&application=NOS.COOPS.TAC.WL&begin_date=$(start_s)&end_date=$(end_s)&datum=MLLW&station=$(station_number)&time_zone=lst_ldt&units=english&format=csv"
|
|
if prediction
|
|
url *= "&interval=hilo"
|
|
end
|
|
tmp_path = download(url)
|
|
df = open(tmp_path) do file
|
|
CSV.read(file, DataFrame)
|
|
end
|
|
dt_local = (dt_s) -> DateTime.(dt_s, "yyyy-mm-dd HH:MM")
|
|
transform!(df,
|
|
"Date Time" => dt_local => "TIME",
|
|
" Water Level" => "WLEVEL")
|
|
|
|
return df[(start_date .<= df.TIME .<= end_date), :]
|
|
end
|
|
|
|
function get_meters_per_pixel(lat, zoom)
|
|
px_per_tile = 256
|
|
return 156543.03 * cos(lat) / (2^zoom)
|
|
end
|
|
|
|
function get_map(lat, lon, zoom)
|
|
px_per_tile = 256
|
|
crop_size = 2 * px_per_tile
|
|
meters_per_pixel = get_meters_per_pixel(lat, zoom)
|
|
xf, yf = get_tile_xy_fraction(lat, lon, zoom)
|
|
x0, y0 = (Int(floor(xf)), Int(floor(yf)))
|
|
xpoint = 1 + px_per_tile + Int(floor((xf - x0) * px_per_tile))
|
|
ypoint = 1 + px_per_tile + Int(floor((yf - y0) * px_per_tile))
|
|
half_size = Int(floor(crop_size / 2))
|
|
ystart = ypoint - half_size
|
|
yend = ypoint + half_size
|
|
xstart = xpoint - half_size
|
|
xend = xpoint + half_size
|
|
xpoint -= xstart
|
|
ypoint -= ystart
|
|
|
|
cache_path = "cache/$(lat)_$(lon)_$zoom.png"
|
|
if isfile(cache_path)
|
|
map = Images.load(cache_path)
|
|
return map, xpoint, ypoint
|
|
end
|
|
|
|
map = Array{RGB{N0f8}, 2}(undef, (px_per_tile * 3, px_per_tile * 3))
|
|
for i in -1:1
|
|
x = x0 + i
|
|
ioffset = 1 + (i + 1) * px_per_tile
|
|
for j in -1:1
|
|
y = y0 + j
|
|
joffset = 1 + (j + 1) * px_per_tile
|
|
img = get_tile(x, y, zoom)
|
|
map[joffset:joffset+px_per_tile-1, ioffset:ioffset+px_per_tile-1] = img
|
|
end
|
|
end
|
|
map = map[ystart:yend, xstart:xend]
|
|
#draw!(map, Ellipse(CirclePointRadius(xpoint, ypoint, 6)),
|
|
# RGB{N0f8}(0.5, 0.0, 0.5))
|
|
Images.save(cache_path, map)
|
|
return map, xpoint, ypoint
|
|
end
|
|
|
|
function get_realtime_data_file(station)
|
|
station = uppercase(station)
|
|
url = "http://www.ndbc.noaa.gov/data/realtime2/$station.txt"
|
|
return download(url)
|
|
end
|
|
|
|
function get_realtime_dataframe(station_name; timezone=tz"America/New_York")
|
|
station_fname = get_realtime_data_file(station_name)
|
|
lines = open(station_fname) do file
|
|
readlines(file)
|
|
end
|
|
dt_array = (y, mon, d, hr, min) -> (astimezone.(ZonedDateTime.(
|
|
DateTime.(y, mon, d, hr, min), tz"UTC"), timezone))
|
|
deg_c_to_f = passmissing(x -> x * 9.0 / 5.0 + 32.0)
|
|
speed_ms_to_mph = passmissing(x -> 0.44704 * x)
|
|
lines[1] = lstrip(lines[1], '#')
|
|
data = join([lines[1:1]; lines[3:end]], "\n")
|
|
df = CSV.read(IOBuffer(data), DataFrame; header=1, delim=" ",
|
|
ignorerepeated=true)
|
|
select!(df,
|
|
[:YY, :MM, :DD, :hh, :mm] => dt_array => :TIME,
|
|
:WDIR => x -> tryparsem.(Int, x),
|
|
:WSPD => x -> speed_ms_to_mph.(tryparsem.(Float32, x)),
|
|
:GST => x -> speed_ms_to_mph.(tryparsem.(Float32, x)),
|
|
:PRES => x -> tryparsem.(Float32, x),
|
|
:ATMP => x -> deg_c_to_f.(tryparsem.(Float32, x)),
|
|
:WTMP => x -> deg_c_to_f.(tryparsem.(Float32, x));
|
|
renamecols=false)
|
|
return df
|
|
end
|
|
|
|
function plot_station_data(station_name, station_map, station_x, station_y,
|
|
df, i)
|
|
fig = Figure(resolution=size(station_map), figure_padding=0)
|
|
ax = CairoMakie.Axis(fig[1, 1])
|
|
hidespines!(ax)
|
|
hidedecorations!(ax)
|
|
plot_station_data_to_axis(ax, station_name, station_map,
|
|
station_x, station_y, df, i)
|
|
return fig
|
|
end
|
|
|
|
function plot_station_data_to_axis(ax, station_name, station_map, station_x, station_y,
|
|
df, i, old_plots=())
|
|
obs_time = df.TIME[i]
|
|
obs_time_s = Dates.format(obs_time, "e, u d H:MM")
|
|
if length(old_plots) == 0
|
|
plot_img = image!(ax, rotr90(station_map))
|
|
else
|
|
plot_img = old_plots[1]
|
|
delete!(ax, old_plots[2])
|
|
delete!(ax, old_plots[3])
|
|
end
|
|
ax.title = "$station_name $obs_time_s"
|
|
wind_rad = passmissing(deg2rad)(90 - df.WDIR[i])
|
|
wind_mph = df.WSPD[i]
|
|
if ismissing(wind_mph) || ismissing(wind_rad)
|
|
wind_mph = 0.0
|
|
wind_rad = 0
|
|
end
|
|
atmp = df.ATMP[i]
|
|
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
|
|
station_y += 50 # Why is this off by 50???
|
|
#water_temp_color = RGB(0, 0, 0.1 + 0.9 * (wtmp - 40) / 40)
|
|
#air_temp_color = RGB(0.5 + (atmp - 40) / 90, 0, 0)
|
|
plot_tooltip = tooltip!(ax, station_x, station_y, txt;
|
|
offset=15, placement=:above)
|
|
plot_arrow = arrows!(ax, [Point2f(station_x, station_y)], [wind_vec],
|
|
fxaa=true, linewidth=2, align=:center, arrowhead='⟰')
|
|
return plot_img, plot_tooltip, plot_arrow
|
|
end
|
|
|
|
function timestamp_range(timestamps, t0::Dates.DateTime)
|
|
return Dates.value.(timestamps) .- Dates.value(t0)
|
|
end
|
|
|
|
struct DateTimeTicks
|
|
t0::Dates.DateTime
|
|
end
|
|
|
|
# See https://github.com/MakieOrg/Makie.jl/issues/442#issuecomment-1446004881
|
|
function Makie.get_ticks(t::DateTimeTicks, any_scale, ::Makie.Automatic, vmin, vmax)
|
|
dateticks, dateticklabels = PlotUtils.optimize_datetime_ticks(
|
|
Dates.value(Dates.DateTime(Dates.Millisecond(Int64(vmin)) + t.t0)),
|
|
Dates.value(Dates.DateTime(Dates.Millisecond(Int64(vmax)) + t.t0)),
|
|
)
|
|
dtformat = "e H:MM"
|
|
labels = [Dates.format(convert(Dates.DateTime, Millisecond(v)), dtformat)
|
|
for v in dateticks]
|
|
return dateticks .- Dates.value(t.t0), labels
|
|
end
|
|
|
|
function save_animation(station_name, out_path, hours, real_hour_per_animation_second)
|
|
info = get_station_info(station_name)
|
|
lat, lon = info.lat, info.lon
|
|
station_map, station_x, station_y = get_map(lat, lon, 14)
|
|
df = get_realtime_dataframe(station_name; timezone=tz"America/New_York")
|
|
diff_hours = (df.TIME[1] - df.TIME[2]).value / 1000 / 3600
|
|
stride = Int(ceil(real_hour_per_animation_second / diff_hours))
|
|
end_index = Int(min(stride * hours + 1, size(df, 1)))
|
|
data_indexes = range(1, end_index; step=stride)
|
|
|
|
df_recent = reverse(@view df[1:end_index, :])
|
|
df_ts = select(df_recent,
|
|
:TIME => x -> Dates.DateTime.(x),
|
|
renamecols=false)
|
|
|
|
df_tides = get_water_levels(info.num, df_ts.TIME[1], df_ts.TIME[end])
|
|
|
|
fig_res = (size(station_map, 1), Int(floor(size(station_map, 2) * 1.5)))
|
|
fig = Figure(resolution=fig_res, figure_padding=0)
|
|
t0 = df_ts.TIME[1]
|
|
ax_temp = CairoMakie.Axis(fig[1, 1], xticks=DateTimeTicks(t0), ylabel="°F")
|
|
ax_wind = CairoMakie.Axis(fig[1, 1], xticks=DateTimeTicks(t0), ylabel="mph",
|
|
yaxisposition=:right, ygridvisible=false,
|
|
xgridvisible=false, xticksvisible=false,
|
|
xticklabelsvisible=false)
|
|
ax_tide= CairoMakie.Axis(fig[2, 1], xticks=DateTimeTicks(t0), ylabel="ft")
|
|
ax = CairoMakie.Axis(fig[3, 1])
|
|
hidespines!(ax)
|
|
hidedecorations!(ax)
|
|
rowsize!(fig.layout, 1, Auto(0.3))
|
|
rowsize!(fig.layout, 2, Auto(0.2))
|
|
|
|
xs_time = timestamp_range(df_ts.TIME, t0)
|
|
lines!(ax_temp, xs_time, df_recent.WTMP, color=:darkblue)
|
|
lines!(ax_temp, xs_time, df_recent.ATMP, color=:brown)
|
|
lines!(ax_wind, xs_time, df_recent.WSPD, color=:skyblue)
|
|
old_vline = vlines!(ax_temp, xs_time[1], color=:black, linewidth=5)
|
|
|
|
lines!(ax_tide, xs_time, df_tides.WLEVEL, color=:blue)
|
|
|
|
old_plots = ()
|
|
record(fig, out_path, data_indexes; framerate=1) do i
|
|
old_plots = plot_station_data_to_axis(ax, station_name, station_map,
|
|
station_x, station_y, df_recent,
|
|
i, old_plots)
|
|
delete!(ax_temp, old_vline)
|
|
old_vline = vlines!(ax_temp, xs_time[i], color=:black, linewidth=5)
|
|
end
|
|
save(replace(out_path, r"\..*" => ".png"), fig)
|
|
end
|
|
|
|
function main()
|
|
if (length(ARGS) < 1)
|
|
println("Usage: create_latest_map.jl station_name [out_path]")
|
|
exit(1)
|
|
end
|
|
station_name = ARGS[1]
|
|
|
|
out_fname = "$station_name.png"
|
|
if length(ARGS) > 1
|
|
out_fname = ARGS[2]
|
|
end
|
|
|
|
if any(endswith.(out_fname, [".gif", ".webm", ".mp4", ".mkv"]))
|
|
save_animation(station_name, out_fname, 48, 1)
|
|
else
|
|
lat, lon = get_station_lat_lon(station_name)
|
|
station_map, station_x, station_y = get_map(lat, lon, 14)
|
|
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
|
|
|
|
main()
|
|
#Profile.print()
|