create_latest_map.jl

@@ -1,291 +1,123 @@
 #!/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")
+	stations = readxml("activestations.xml")
-    station_node = findfirst("/stations/station[@id=\"$station_id\"]", stations)
+	station_node = findfirst("/stations/station[@id=\"$station_id\"]", stations)
-    return parse.(Float32, [station_node["lat"], station_node["lon"]])
+	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
+	n = 2^zoom
-    lat_rad = lat * pi / 180
+	lat_rad = lat * pi / 180
-    xtile = (lon + 180) / 360 * n
+	xtile = (lon + 180) / 360 * n
-        ytile = (1 - log(tan(lat_rad) + 1 / cos(lat_rad)) / pi) / 2 * n
+    ytile = (1 - log(tan(lat_rad) + 1 / cos(lat_rad)) / pi) / 2 * n
-    return xtile, ytile
+	return xtile, ytile
 end
 
 function get_tile_xy(lat, lon, zoom)
-    xf, yf = get_tile_xy_fraction(lat, lon, zoom)
+	xf, yf = get_tile_xy_fraction(lat, lon, zoom)
-    return (Int(floor(xf)), Int(floor(yf)))
+	return (Int(floor(xf)), Int(floor(yf)))
 end
 
 function get_tile(x, y, zoom)
-    url = "https://tile.openstreetmap.org/$zoom/$x/$y.png"
+	url = "https://tile.openstreetmap.org/$zoom/$x/$y.png"
-    tmp_path = download(url)
+	tmp_path = download(url)
-    tile = load(tmp_path)
+	tile = load(tmp_path)
-    return tile
+	return tile
 end
 
 function get_tile_by_lat_lon(lat, lon, zoom)
-    xtile, ytile = get_tile_xy(lat, lon, zoom)
+	xtile, ytile = get_tile_xy(lat, lon, zoom)
-    return get_tile(xtile, ytile, 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
+	px_per_tile = 256
-    return 156543.03 * cos(lat) / (2^zoom)
+	return 156543.03 * cos(lat) / (2^zoom)
 end
 
 function get_map(lat, lon, zoom)
-    px_per_tile = 256
+	px_per_tile = 256
-    crop_size = 2 * px_per_tile
+	meters_per_pixel = get_meters_per_pixel(lat, zoom)
-    meters_per_pixel = get_meters_per_pixel(lat, zoom)
+	xf, yf = get_tile_xy_fraction(lat, lon, zoom)
-    xf, yf = get_tile_xy_fraction(lat, lon, zoom)
+	x0, y0 = (Int(floor(xf)), Int(floor(yf)))
-    x0, y0 = (Int(floor(xf)), Int(floor(yf)))
+	xpoint = 1 + px_per_tile + Int(floor((xf - x0) * px_per_tile))
-    xpoint = 1 + px_per_tile + Int(floor((xf - x0) * px_per_tile))
+	ypoint = 1 + px_per_tile + Int(floor((yf - y0) * px_per_tile))
-    ypoint = 1 + px_per_tile + Int(floor((yf - y0) * px_per_tile))
+	map = Array{RGB{N0f8}, 2}(undef, (px_per_tile * 3, px_per_tile * 3))
-    half_size = Int(floor(crop_size / 2))
+	for i in -1:1
-    ystart = ypoint - half_size
+		x = x0 + i
-    yend = ypoint + half_size
+		ioffset = 1 + (i + 1) * px_per_tile
-    xstart = xpoint - half_size
+		for j in -1:1
-    xend = xpoint + half_size
+			y = y0 + j
-    xpoint -= xstart
+		    joffset = 1 + (j + 1) * px_per_tile
-    ypoint -= ystart
+			img = get_tile(x, y, zoom)
-
+			map[joffset:joffset+px_per_tile-1, ioffset:ioffset+px_per_tile-1] = img
-    cache_path = "cache/$(lat)_$(lon)_$zoom.png"
+		end
-    if isfile(cache_path)
+	end
-        map = Images.load(cache_path)
+	draw!(map, Ellipse(CirclePointRadius(xpoint, ypoint, 6)),
-        return map, xpoint, ypoint
+		  RGB{N0f8}(0.5, 0.0, 0.5))
-    end
+	return map, xpoint, ypoint
-
-    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"
+	url = "http://www.ndbc.noaa.gov/data/realtime2/$station.txt"
-    return download(url)
+	return download(url)
 end
 
 function get_realtime_dataframe(station_name; timezone=tz"America/New_York")
-    station_fname = get_realtime_data_file(station_name)
+	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)
+	deg_c_to_f = passmissing(x -> x * 9.0 / 5.0 + 32.0)
-    speed_ms_to_mph = passmissing(x -> 0.44704 * x)
+	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,
+	select!(df,
-            [:YY, :MM, :DD, :hh, :mm] => dt_array => :TIME,
+		[:YY, :MM, :DD, :hh, :mm] => dt_array => :TIME,
-            :WDIR => x -> tryparsem.(Int, x),
+		:WDIR => x -> tryparsem.(Int, x),
-            :WSPD => x -> speed_ms_to_mph.(tryparsem.(Float32, x)),
+		:WSPD => x -> speed_ms_to_mph.(tryparsem.(Float32, x)),
-            :GST => x -> speed_ms_to_mph.(tryparsem.(Float32, x)),
+		:GST => x -> speed_ms_to_mph.(tryparsem.(Float32, x)),
-            :PRES => x -> tryparsem.(Float32, x),
+		:PRES => x -> tryparsem.(Float32, x),
-            :ATMP => x -> deg_c_to_f.(tryparsem.(Float32, x)),
+		:ATMP => x -> deg_c_to_f.(tryparsem.(Float32, x)),
-            :WTMP => x -> deg_c_to_f.(tryparsem.(Float32, x));
+		:WTMP => x -> deg_c_to_f.(tryparsem.(Float32, x));
-            renamecols=false)
+		renamecols=false)
-    return df
+	return df
-end
+end
-
+
-function plot_station_data(station_name, station_map, station_x, station_y,
+function plot_station_data(station_name, station_map, station_x, station_y, df, i)
-                           df, i)
+	obs_time = df.TIME[i]
-    fig = Figure(resolution=size(station_map), figure_padding=0)
+	fig, ax, p = image(rotr90(station_map),
-    ax = CairoMakie.Axis(fig[1, 1])
+		               axis=(aspect=DataAspect(),
-    hidespines!(ax)
+                             title="$station_name $obs_time"))
-    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))
+	wind_rad = deg2rad(90 - df.WDIR[1])
-    rowsize!(fig.layout, 2, Auto(0.2))
+	wind_mph = df.WSPD[i]
-
+	txt = @sprintf "Air: %.1f F, Water: %.1f F, Wind: %.1f mph" df.ATMP[i] df.WTMP[i] wind_mph
-    xs_time = timestamp_range(df_ts.TIME, t0)
+	wind_vec = 50 * Vec2f(cos(wind_rad), sin(wind_rad)) * wind_mph / 10
-    lines!(ax_temp, xs_time, df_recent.WTMP, color=:darkblue)
+	tooltip!(ax, station_x, station_y, txt; offset=10)
-    lines!(ax_temp, xs_time, df_recent.ATMP, color=:brown)
+	arrows!(ax, [Point2f(station_x, station_y)], [wind_vec], fxaa=true,
-    lines!(ax_wind, xs_time, df_recent.WSPD, color=:skyblue)
+	        linewidth=2, align=:center, arrowhead='⟰')
-    old_vline = vlines!(ax_temp, xs_time[1], color=:black, linewidth=5)
+	return fig
-
-    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()
@@ -296,21 +128,16 @@ 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"]))
+    lat, lon = get_station_lat_lon(station_name)
-        save_animation(station_name, out_fname, 48, 1)
+	station_map, station_x, station_y = get_map(lat, lon, 14)
-    else
+    df = get_realtime_dataframe(station_name; timezone=tz"America/New_York")
-        lat, lon = get_station_lat_lon(station_name)
+    fig = plot_station_data(station_name, station_map, station_x, station_y,
-        station_map, station_x, station_y = get_map(lat, lon, 14)
+                            df, 1)
-        df = get_realtime_dataframe(station_name; timezone=tz"America/New_York")
+    save("$station_name.png", fig)
-        fig = plot_station_data(station_name, station_map, station_x, station_y,
-                                df, 1)
-        save(out_fname, fig)
-    end
 end
 
 main()
-#Profile.print()