julia: add tide data

LICENSE

@@ -1,23 +1,26 @@
-Copyright <yyyy, yyyy> The Open Group
+Copyright 2023 Bryce Allen <bdallen@uchicago.edu>
 
-Permission to use, copy, modify, distribute, and sell this software and
-its documentation for any purpose is hereby granted without fee,
-provided that the above copyright notice appear in all copies and that
-both that copyright notice and this permission notice appear in
-supporting documentation.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
 
-The above copyright notice and this permission notice shall be included
-in all copies or substantial portions of the Software.
+1. Redistributions of source code must retain the above copyright notice, this
+list of conditions and the following disclaimer.
 
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
-OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-IN NO EVENT SHALL BE LIABLE FOR ANY CLAIM, DAMAGES
-OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
-OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
-THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+2. Redistributions in binary form must reproduce the above copyright notice,
+this list of conditions and the following disclaimer in the documentation
+and/or other materials provided with the distribution.
 
-Except as contained in this notice, the name of The Open Group
-shall not be used in advertising or otherwise to promote the sale, use
-or other dealings in this Software without prior written authorization
-from The Open Group.
+3. Neither the name of the copyright holder nor the names of its contributors
+may be used to endorse or promote products derived from this software without
+specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

create_latest_map.jl

@@ -1,15 +1,38 @@
 #!/usr/bin/env julia
-#
+
 using DataFrames
 using Dates
 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
 
@@ -103,23 +168,126 @@ function get_realtime_dataframe(station_name; timezone=tz"America/New_York")
     return df
 end
 
-function plot_station_data(station_name, station_map, station_x, station_y, df, i)
-	obs_time = df.TIME[i]
-	fig, ax, p = image(rotr90(station_map),
-		               axis=(aspect=DataAspect(),
-                             title="$station_name $obs_time"))
+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)
-	wind_rad = deg2rad(90 - df.WDIR[1])
-	wind_mph = df.WSPD[i]
-	txt = @sprintf "Air: %.1f F, Water: %.1f F, Wind: %.1f mph" df.ATMP[i] df.WTMP[i] 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='⟰')
+    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]")
@@ -128,16 +296,21 @@ function main()
     station_name = ARGS[1]
 
     out_fname = "$station_name.png"
-    if (length(ARGS) > 1)
+    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("$station_name.png", fig)
+        save(out_fname, fig)
+    end
 end
 
 main()
+#Profile.print()