diff --git a/create_latest_map.jl b/create_latest_map.jl new file mode 100644 index 0000000..575982d --- /dev/null +++ b/create_latest_map.jl @@ -0,0 +1,143 @@ +#!/usr/bin/env julia +# +using DataFrames +using Dates +using CSV +using EzXML +using Images +using Colors +using ImageDraw +using CairoMakie +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"]]) +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_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 + 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) +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) + obs_time = df.TIME[i] + fig, ax, p = image(rotr90(station_map), + axis=(aspect=DataAspect(), + title="$station_name $obs_time")) + 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='⟰') + return 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 + + 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) +end + +main() diff --git a/notebook.jl b/notebook.jl index 0f4595f..d265376 100644 --- a/notebook.jl +++ b/notebook.jl @@ -141,13 +141,16 @@ function plot_station_data(station_map, station_x, station_y, df, i) obs_time = df.TIME[i] fig, ax, p = image(rotr90(station_map), axis=(aspect=DataAspect(), title="DUKN7 $obs_time")) + hidespines!(ax) + hidedecorations!(ax) wind_rad = deg2rad(90 - df.WDIR[1]) - wind_mph = df.WSPD[i] + wind_mph = 0.1#df.WSPD[i] txt = @sprintf "Air: %.1f F, Water: %.1f F, Wind: %.1f mph" df.ATMP[i] df.WTMP[i] wind_mph - wind_vec = 20 * Vec2f(cos(wind_rad), sin(wind_rad)) * wind_mph / 10 - tooltip!(station_x, station_y, txt) - arrows!([Point2f(station_x, station_y)], [wind_vec], fxaa=true, - linewidth=2, align=:center) + 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, lengthscale=1, arrowhead='⟰') + # ⋀ fig end