diff --git a/create_latest_map.jl b/create_latest_map.jl index 5439dda..e8117b7 100644 --- a/create_latest_map.jl +++ b/create_latest_map.jl @@ -1,5 +1,5 @@ #!/usr/bin/env julia -# + using DataFrames using Dates using CSV @@ -12,140 +12,157 @@ using Printf using TimeZones 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"]]) + 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 + 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))) + 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 + 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) + xtile, ytile = get_tile_xy(lat, lon, zoom) + return get_tile(xtile, ytile, zoom) end function get_meters_per_pixel(lat, zoom) - px_per_tile = 256 - return 156543.03 * cos(lat) / (2^zoom) + px_per_tile = 256 + return 156543.03 * cos(lat) / (2^zoom) end function get_map(lat, lon, zoom) - px_per_tile = 256 - 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)) - 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 - draw!(map, Ellipse(CirclePointRadius(xpoint, ypoint, 6)), - RGB{N0f8}(0.5, 0.0, 0.5)) - return map, xpoint, ypoint + px_per_tile = 256 + 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)) + 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 + draw!(map, Ellipse(CirclePointRadius(xpoint, ypoint, 6)), + RGB{N0f8}(0.5, 0.0, 0.5)) + 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) + 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) + 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) + 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 + 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)) - 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 + fig = Figure(resolution=size(station_map)) + 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) - obs_time = df.TIME[i] - obs_time_s = Dates.format(obs_time, "e, u d H:MM") - image!(ax, rotr90(station_map)) - ax.title = "$station_name $obs_time_s" - 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='⟰') - return ax + obs_time = df.TIME[i] + obs_time_s = Dates.format(obs_time, "e, u d H:MM") + image!(ax, rotr90(station_map)) + ax.title = "$station_name $obs_time_s" + wind_rad = deg2rad(90 - df.WDIR[1]) + 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 + 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 end -function save_animation(station_name, out_path, days, real_hour_per_animation_second=1) +function save_animation(station_name, out_path, hours, real_hour_per_animation_second) 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") - diff_hours = (df.TIME[0] - df.TIME[1]).value / 1000 / 3600 - fps = real_hour_per_animation_second / diff_hours - nframes = fps * days * 24 * 3600 + 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) fig = Figure(resolution=size(station_map)) - ax = Axis(fig[1, 1]) + ax = CairoMakie.Axis(fig[1, 1]) hidespines!(ax) hidedecorations!(ax) - record(fig, out_path, range(1, nframes); framerate=fps) do i + 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, nframes - i + 1) + station_x, station_y, df, end_index - i + 1) end end @@ -157,16 +174,20 @@ function main() station_name = ARGS[1] out_fname = "$station_name.png" - if (length(ARGS) > 1) + if length(ARGS) > 1 out_fname = ARGS[2] end - 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) + if any(endswith.(out_fname, [".gif", ".webm", ".mp4", ".mkv"])) + save_animation(station_name, out_fname, 24, 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()