bda
/
plot-buoy-condition
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

julia: add tide data

Browse Source
main
Bryce Allen 2 years ago
parent 7173dea822
commit 7e63c8901b
1 changed files with 145 additions and 22 deletions
Show Stats Download Patch File Download Diff File

167
create_latest_map.jl
Unescape View File

@ -6,10 +6,33 @@ using CSV
using EzXML
using Images
using Colors
using ImageDraw
#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")
@ -44,6 +67,31 @@ function get_tile_by_lat_lon(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)
@ -51,11 +99,26 @@ 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
@ -67,8 +130,10 @@ function get_map(lat, lon, zoom)
map[joffset:joffset+px_per_tile-1, ioffset:ioffset+px_per_tile-1] = img
end
end
draw!(map, Ellipse(CirclePointRadius(xpoint, ypoint, 6)),
RGB{N0f8}(0.5, 0.0, 0.5))
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
@ -105,7 +170,7 @@ end
function plot_station_data(station_name, station_map, station_x, station_y,
df, i)
fig = Figure(resolution=size(station_map))
fig = Figure(resolution=size(station_map), figure_padding=0)
ax = CairoMakie.Axis(fig[1, 1])
hidespines!(ax)
hidedecorations!(ax)
@ -115,12 +180,18 @@ function plot_station_data(station_name, station_map, station_x, station_y,
end
function plot_station_data_to_axis(ax, station_name, station_map, station_x, station_y,
df, i)
df, i, old_plots=())
obs_time = df.TIME[i]
obs_time_s = Dates.format(obs_time, "e, u d H:MM")
image!(ax, rotr90(station_map))
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 = deg2rad(90 - df.WDIR[1])
wind_rad = passmissing(deg2rad)(90 - df.WDIR[i])
wind_mph = df.WSPD[i]
if ismissing(wind_mph) || ismissing(wind_rad)
wind_mph = 0.0
@ -136,34 +207,85 @@ function plot_station_data_to_axis(ax, station_name, station_map, station_x, sta
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
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)
lat, lon = get_station_lat_lon(station_name)
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))
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)
data_indexes = range(1, end_index; step=stride)
fig = Figure(resolution=size(station_map))
ax = CairoMakie.Axis(fig[1, 1])
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
plot_station_data_to_axis(ax, station_name, station_map,
station_x, station_y, df, end_index - i + 1)
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()
@ -179,7 +301,7 @@ function main()
end
if any(endswith.(out_fname, [".gif", ".webm", ".mp4", ".mkv"]))
save_animation(station_name, out_fname, 24, 1)
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)
@ -191,3 +313,4 @@ function main()
end
main()
#Profile.print()

Write Preview
Loading…
Cancel
Save