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

@@ -0,0 +1,316 @@
+#!/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()