diff --git a/create_latest_map.jl b/create_latest_map.jl index e8117b7..6f952d7 100644 --- a/create_latest_map.jl +++ b/create_latest_map.jl @@ -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) + + # + 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()