#!/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()