julia: hide axis, add script

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()

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
+	wind_vec = 50 * Vec2f(cos(wind_rad), sin(wind_rad)) * wind_mph / 10
-	tooltip!(station_x, station_y, txt)
+	tooltip!(ax, station_x, station_y, txt; offset=10)
-	arrows!([Point2f(station_x, station_y)], [wind_vec], fxaa=true,
+	arrows!(ax, [Point2f(station_x, station_y)], [wind_vec], fxaa=true,
-	        linewidth=2, align=:center)
+	        linewidth=2, align=:center, lengthscale=1, arrowhead='⟰')
+	# ⋀
 	fig
 end