leaflet_svelte/docs/wind-vis-math.tex
2026-06-17 00:20:55 +09:00

395 lines
14 KiB
TeX

\documentclass[11pt,a4paper]{article}
\usepackage{siunitx}
\usepackage{amsmath,amssymb,amsthm}
\usepackage{geometry}
\usepackage{hyperref}
\usepackage{booktabs}
\geometry{margin=2.5cm}
\title{Wind Visualisation -- Mathematical Reference\\
\large stratoflights-predictor frontend}
\author{stratoflights}
\date{}
\begin{document}
\maketitle
\tableofcontents
\bigskip
%---------------------------------------------------------------------------
\section{Wind Field Grid Format}
%---------------------------------------------------------------------------
The predictor's \texttt{GET /api/v1/wind/field} endpoint returns a
two-element JSON array $[C_U,\,C_V]$, each element having the shape
\begin{verbatim}
{ "header": { "nx", "ny", "lo1", "la1", "lo2", "la2",
"dx", "dy", "refTime", ... },
"data": [ float, ... ] // flat row-major array, length = nx * ny
}
\end{verbatim}
\noindent where $C_U$ contains the \emph{eastward} (zonal) component $u$
and $C_V$ contains the \emph{northward} (meridional) component $v$, both
in \si{m\,s^{-1}}.
\subsection{Grid coordinates}
Let $n_x$ and $n_y$ be the number of grid points along the longitude and
latitude axes respectively. The coordinate of grid cell $(i,\,j)$
($i = 0, \ldots, n_x-1$;\; $j = 0, \ldots, n_y-1$) is
\begin{align}
\lambda_{i} &= \operatorname{wrap}\!\left(\lambda_1 + i\,\Delta\lambda\right), \label{eq:lng}\\
\varphi_{j} &= \varphi_1 + j\,\Delta\varphi, \label{eq:lat}
\end{align}
\noindent where $\lambda_1 = \texttt{lo1}$ and $\varphi_1 = \texttt{la1}$.
\paragraph{Increments from the extent, not \texttt{dx}/\texttt{dy}.}
The predictor reports \texttt{dx} and \texttt{dy} as positive magnitudes
\emph{regardless of scan direction}, and emits longitudes in the
$[0,360)$ convention. The per-step increments are therefore derived from
the grid extent so the last row/column lands exactly on
$(\texttt{la2},\texttt{lo2})$ and the scan direction follows automatically:
\begin{equation}
\Delta\lambda = \frac{(\texttt{lo2}-\texttt{lo1}) \bmod 360}{n_x-1},
\qquad
\Delta\varphi = \frac{\texttt{la2}-\texttt{la1}}{n_y-1}.
\label{eq:increments}
\end{equation}
For the standard GFS global grid $\varphi_1 = 90^\circ$, $\varphi_2 = -90^\circ$,
so $\Delta\varphi = -10^\circ$ (rows scan southward) without any case
distinction. (When $n_x=1$ or $n_y=1$ the magnitude \texttt{dx}/\texttt{dy}
is used as a fallback.)
\paragraph{Longitude wrapping.}
Because longitudes are emitted in $[0,360)$ but the map renders in
$(-180,180]$, each longitude is wrapped:
\begin{equation}
\operatorname{wrap}(\lambda) = \big((\lambda + 180) \bmod 360\big) - 180.
\label{eq:wrap}
\end{equation}
Without \eqref{eq:wrap} a grid sampled near the prime meridian
(e.g.\ $\texttt{lo1}=358$ for a query at $-2^\circ$) would be rendered a
full $360^\circ$ away from the trajectory.
The flat data index for cell $(i,\,j)$ is
\begin{equation}
k = j\,n_x + i.
\label{eq:index}
\end{equation}
%---------------------------------------------------------------------------
\section{Bilinear Interpolation of Wind Components}
%---------------------------------------------------------------------------
The particle renderer samples the wind at arbitrary points (the unprojected
pixel positions of individual particles), so bilinear interpolation of the
$[u,v]$ grid is performed at run time for every particle, every frame.
Given a query point $(\lambda,\,\varphi)$ inside the grid, locate the
surrounding four cells
\begin{align*}
i_0 &= \left\lfloor \frac{\lambda - \lambda_1}{\Delta\lambda} \right\rfloor,
&
j_0 &= \left\lfloor \frac{\varphi - \varphi_1}{\Delta\varphi} \right\rfloor,
\end{align*}
and define the fractional offsets
\begin{equation*}
s = \frac{\lambda - \lambda_{i_0}}{\Delta\lambda},
\qquad
t = \frac{\varphi - \varphi_{j_0}}{\Delta\varphi},
\qquad
s,t \in [0,1].
\end{equation*}
The bilinearly interpolated value of any scalar field $f$ at $(\lambda,\varphi)$
is
\begin{equation}
f(\lambda,\varphi)
= (1-s)(1-t)\,f_{i_0,j_0}
+ s(1-t)\,f_{i_0+1,j_0}
+ (1-s)t\,f_{i_0,j_0+1}
+ st\,f_{i_0+1,j_0+1}.
\label{eq:bilinear}
\end{equation}
Applied independently to $u$ and $v$, equation~\eqref{eq:bilinear} gives
the interpolated wind vector at any interior point.
%---------------------------------------------------------------------------
\section{Particle Advection and Rendering}
%---------------------------------------------------------------------------
The wind is visualised as a dense field of particles that flow with the wind,
in the style of \texttt{leaflet-velocity} / cambecc's \emph{earth}. Particles
live in screen (CSS-pixel) space; each animation frame they are advected by
the local wind and drawn as short fading trails.
\subsection{Speed (magnitude)}
The wind speed used for colouring is the Euclidean norm of the horizontal
wind vector returned by the interpolator~\eqref{eq:bilinear}:
\begin{equation}
\lVert\mathbf{w}\rVert = \sqrt{u^2 + v^2}.
\label{eq:speed}
\end{equation}
\subsection{Advection through the map projection}
A particle at pixel position $\mathbf{x} = (x,y)$ is unprojected to geographic
coordinates $(\lambda,\varphi) = P^{-1}(\mathbf{x})$, where $P$ is the
MapLibre projection (Web Mercator composed with the current view transform).
The wind $\mathbf{w}=(u,v)$ in the east--north frame must be expressed in
pixel space; this is done with the local Jacobian of $P$, estimated by finite
differences with a small $\varepsilon$ (degrees):
\begin{align}
\mathbf{J}_\lambda &= \frac{P(\lambda+\varepsilon,\varphi) - \mathbf{x}}{\varepsilon},
&
\mathbf{J}_\varphi &= \frac{P(\lambda,\varphi+\varepsilon) - \mathbf{x}}{\varepsilon}.
\label{eq:jacobian}
\end{align}
\noindent The pixel-space velocity is the Jacobian applied to the wind vector,
scaled by a dimensionless speed factor $c$ (the configurable
\texttt{particleSpeed} times a base constant):
\begin{equation}
\dot{\mathbf{x}} = c\,\big(\mathbf{J}_\lambda\, u + \mathbf{J}_\varphi\, v\big).
\label{eq:pixel_velocity}
\end{equation}
Because $\mathbf{J}_\varphi$ points toward $-y$ (north is up), a northward wind
moves the particle up the screen, and the projection automatically supplies
the correct scale at every zoom level and the meridian-convergence distortion
near the poles. The particle is integrated one explicit Euler step per frame:
\begin{equation}
\mathbf{x}_{t+1} = \mathbf{x}_t + \dot{\mathbf{x}}.
\label{eq:euler}
\end{equation}
A particle whose position leaves the data grid (interpolator returns
\texttt{null}) or whose age exceeds $A_{\max}$ frames is respawned at a random
pixel that has wind, with a randomised initial age to desynchronise the
population.
\subsection{Colour mapping}
Speed is mapped to a 15-stop colour ramp (blue $\to$ red) by linear
quantisation into the range $[v_{\min}, v_{\max}]$:
\begin{equation}
k(\lVert\mathbf{w}\rVert) =
\operatorname{round}\!\left(
\frac{\lVert\mathbf{w}\rVert - v_{\min}}{v_{\max}-v_{\min}}\,(K-1)
\right),
\quad k \in \{0,\dots,K-1\},
\end{equation}
\noindent clamped to the ramp, where $K=15$ and $v_{\max}$ is the configurable
\texttt{maxVelocity}. Trails sharing a colour bucket are batched into a single
stroked path per frame.
\subsection{Trail fading}
Rather than clearing the canvas each frame, the previous frame is faded toward
\emph{transparency} (so the basemap stays visible) by compositing a translucent
rectangle with the \texttt{destination-in} operator:
\begin{equation}
\alpha_{t+1}(\mathbf{x}) = \rho\,\alpha_t(\mathbf{x}),
\qquad \rho \in [0,1),
\end{equation}
\noindent where $\rho$ is the configurable \texttt{trailPersistence}; an
existing trail therefore decays geometrically, retaining a visible tail of
length $\sim 1/(1-\rho)$ frames. New trail segments are then drawn with the
\texttt{source-over} operator.
\subsection{Particle count}
The particle population is proportional to the canvas area:
\begin{equation}
N = \min\!\big(N_{\max},\; W H \, \mu \, \texttt{density}\big),
\end{equation}
\noindent with base multiplier $\mu = 1/350$ particles per pixel, the
configurable \texttt{density} factor, and a hard cap $N_{\max}=6000$ to bound
per-frame cost. The field is re-seeded on resize and on map \texttt{moveend};
trails are cleared while the map is panning/zooming and resume afterwards.
%---------------------------------------------------------------------------
\section{Trajectory Bounding Box}
%---------------------------------------------------------------------------
Let the prediction trajectory be the sequence of points
$\{(\varphi_k, \lambda_k)\}_{k=0}^{N-1}$. The axis-aligned bounding box
with margin $m$ (degrees) is
\begin{align}
\varphi_{\min}' &= \min_k \varphi_k - m, &
\varphi_{\max}' &= \max_k \varphi_k + m, \notag\\
\lambda_{\min}' &= \min_k \lambda_k - m, &
\lambda_{\max}' &= \max_k \lambda_k + m.
\label{eq:bbox}
\end{align}
A wind field request is suppressed when either span exceeds the configured
maximum $D_{\max}$:
\begin{equation}
(\varphi_{\max}' - \varphi_{\min}') > D_{\max}
\quad\text{or}\quad
(\lambda_{\max}' - \lambda_{\min}') > D_{\max}.
\label{eq:bbox_guard}
\end{equation}
Default values: $m = 1^\circ$, $D_{\max} = 20^\circ$.
%---------------------------------------------------------------------------
\section{Trajectory Altitude Lookup}
%---------------------------------------------------------------------------
The prediction result contains two parallel arrays:
$\mathbf{p} = \{(\varphi_k, \lambda_k, a_k)\}$ (flight path with altitude
in metres) and $\mathbf{t} = \{t_k\}$ (absolute epoch timestamps in
milliseconds, $t_k < t_{k+1}$).
Given a flight-time offset $\delta$ (milliseconds from launch), the
absolute query time is $T = t_0 + \delta$. The index of the immediately
following trajectory point is found by binary search:
\begin{equation}
k^* = \min\{k \in \{0,\ldots,N-1\} : t_k \ge T\}.
\label{eq:binsearch}
\end{equation}
The trajectory altitude used for the wind query is then
$a_{k^*}$ (nearest-neighbour in time, no interpolation needed since the
pre-fetch frames are already spaced far apart relative to the step size).
%---------------------------------------------------------------------------
\section{Time Series Pre-Fetching}
%---------------------------------------------------------------------------
\subsection{Frame schedule}
Let $F$ be the total flight duration in milliseconds and $\Delta T$ the
pre-fetch interval (milliseconds, default $15 \times 60\,000$). The set
of pre-fetch time offsets is
\begin{equation}
\mathcal{S} = \{0,\;\Delta T,\;2\Delta T,\;\ldots,\;F\},
\label{eq:schedule}
\end{equation}
where the last element is clamped to $F$ (so the landing point is always
included). For each $\delta \in \mathcal{S}$ a wind field is fetched at
absolute time $T = t_0 + \delta$ and trajectory altitude $a_{k^*(\delta)}$
within the bounding box~\eqref{eq:bbox}.
\subsection{Guard conditions}
Pre-fetching is skipped entirely when
\begin{enumerate}
\item $F > F_{\max}$ (flight too long; default $F_{\max} = 4$\,h), or
\item either bounding-box dimension exceeds $D_{\max}$
(equation~\eqref{eq:bbox_guard}).
\end{enumerate}
Both limits are configurable in the application settings.
\subsection{Temporal interpolation between frames}
During playback at flight-time offset $\delta$, rather than snapping to the
nearest cached frame (which would make the wind jump every $\Delta T$), the
displayed field is \emph{linearly interpolated} between the two bracketing
frames $\delta_p \le \delta < \delta_{p+1}$ ($\delta_p,\delta_{p+1}\in\mathcal{S}$).
With blend factor
\begin{equation}
\alpha = \frac{\delta - \delta_p}{\delta_{p+1} - \delta_p} \in [0,1),
\label{eq:blend_alpha}
\end{equation}
each wind component is blended cell-for-cell:
\begin{equation}
u(\delta) = (1-\alpha)\,u^{(p)} + \alpha\,u^{(p+1)},
\qquad
v(\delta) = (1-\alpha)\,v^{(p)} + \alpha\,v^{(p+1)}.
\label{eq:time_lerp}
\end{equation}
This is valid because every frame is fetched over the \emph{same} bounding
box~\eqref{eq:bbox} and step, so the grids are cell-aligned ($n_x,n_y,
\lambda_1,\varphi_1$ identical) and~\eqref{eq:time_lerp} requires only an
element-wise blend of the two flat $[u,v]$ arrays. Note that the two frames
are sampled at slightly different altitudes $a_{k^*(\delta_p)}$ and
$a_{k^*(\delta_{p+1})}$; the blend therefore also interpolates across the
balloon's changing altitude, which is the desired behaviour. The result is
continuous wind evolution along the route while still requiring only
$|\mathcal{S}|$ network requests per trajectory.
The blended field~\eqref{eq:time_lerp} is handed to the particle renderer
(\S\,Particle Advection) as the interpolator swapped in each frame, so the
flowing particles advect through the smoothly time-evolving wind without any
visible stepping between pre-fetch frames.
%---------------------------------------------------------------------------
\section{Complexity and Performance Notes}
%---------------------------------------------------------------------------
\subsection{Grid size and render cost}
For a regional bounding box of $L_\varphi \times L_\lambda$ degrees and
step $h$ degrees, the number of grid cells fetched is
\begin{equation}
n = \left\lceil \frac{L_\varphi}{h} \right\rceil
\left\lceil \frac{L_\lambda}{h} \right\rceil.
\end{equation}
With the default trajectory step $h = 1^\circ$ and worst-case dimensions
$20^\circ \times 20^\circ$ this is $n = 400$ cells; the static global view at
$h = 2^\circ$ is $90 \times 180 = 16\,200$ cells. Crucially, the grid size
only affects the \emph{fetch} and the per-particle bilinear lookup, not the
render budget: the per-frame cost is governed by the particle count $N$
(eq.~for $N$), capped at $N_{\max}=6000$, independent of the grid resolution
or the geographic extent. Each particle costs one unprojection, one
interpolation, and two projections (for the Jacobian) per frame.
\subsection{Request count}
The total number of requests for one trajectory is
\begin{equation}
|\mathcal{S}| = \left\lfloor \frac{F}{\Delta T} \right\rfloor + 1
\;\le\; \frac{F_{\max}}{\Delta T} + 1.
\end{equation}
With the defaults $F_{\max} = 4$\,h and $\Delta T = 15$\,min:
$|\mathcal{S}| \le 17$. Requests are issued sequentially to avoid
bursty load on the predictor.
\subsection{Caching}
The in-memory cache (keyed by the full parameter tuple) ensures that
re-running a prediction with the same parameters, scrubbing the timeline
back and forth, or toggling the wind layer all reuse existing responses.
\end{document}