ac9f9da3c4
map_h was subtracting MARGIN twice but the map_area only has a single bottom margin, leaving the last 8 rows of the land background unrendered. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
337 lines
12 KiB
Rust
337 lines
12 KiB
Rust
use crate::projection::OrthoProjection;
|
||
use anyhow::{Context, Result};
|
||
use font8x8::{UnicodeFonts, BASIC_FONTS};
|
||
use plotters::prelude::*;
|
||
use rayon::prelude::*;
|
||
use spade::{DelaunayTriangulation, HasPosition, Point2, Triangulation};
|
||
use std::path::Path;
|
||
|
||
const IMG_WIDTH: u32 = 900;
|
||
const IMG_HEIGHT: u32 = 600;
|
||
|
||
/// Map image centre and view window (must match what the Python script used).
|
||
const CENTER_LAT: f64 = 52.56;
|
||
const CENTER_LON: f64 = 13.08;
|
||
const HALF_W_DEG: f64 = 0.8; // horizontal half-extent in degrees of arc
|
||
const HALF_H_DEG: f64 = 0.4; // vertical half-extent
|
||
|
||
const BACKGROUND: RGBColor = RGBColor(240, 248, 255); // pale sky-blue
|
||
const LAND_BG: RGBColor = RGBColor(142, 178, 35); // strong green for land
|
||
const BLACK: RGBColor = RGBColor( 0, 0, 0);
|
||
const RED: RGBColor = RGBColor(220, 30, 30);
|
||
|
||
/// Locations to mark on the map: (label, longitude, latitude)
|
||
const LOCATIONS: &[(&str, f64, f64)] = &[
|
||
("Falkensee", 13.08, 52.56),
|
||
("Pausin", 13.04, 52.64),
|
||
("Nauen", 12.88, 52.605),
|
||
("Hennigsdorf", 13.205, 52.64),
|
||
("Ketzin/Brueckenkopf", 12.82, 52.49),
|
||
("Potsdam", 13.06, 52.40),
|
||
("Berlin (Mitte)", 13.39, 52.52),
|
||
];
|
||
|
||
#[derive(Clone, Copy)]
|
||
struct CloudVertex {
|
||
position: Point2<f64>,
|
||
cloud: f32,
|
||
}
|
||
|
||
impl HasPosition for CloudVertex {
|
||
type Scalar = f64;
|
||
fn position(&self) -> Point2<f64> {
|
||
self.position
|
||
}
|
||
}
|
||
|
||
/// Map cloud cover (0–100 %) to a colour by blending the land background toward
|
||
/// white. Returns `None` for near-zero cover so the land background shows through.
|
||
fn cloud_color(cover: f32) -> Option<RGBColor> {
|
||
if cover < 1.0 {
|
||
return None;
|
||
}
|
||
let t = (cover / 100.0).clamp(0.0, 1.0);
|
||
let blend = |base: u8| -> u8 { (base as f32 + t * (255.0 - base as f32)).round() as u8 };
|
||
Some(RGBColor(blend(LAND_BG.0), blend(LAND_BG.1), blend(LAND_BG.2)))
|
||
}
|
||
|
||
/// Draw ASCII text using the built-in 8×8 pixel font.
|
||
///
|
||
/// Coordinates are in the drawing area's local coordinate system.
|
||
/// `scale` is the pixel-multiplier: 1 → 8×8 px/glyph, 2 → 16×16, etc.
|
||
/// Characters outside Basic Latin are rendered as '?'.
|
||
fn draw_pixel_text(
|
||
area: &DrawingArea<BitMapBackend, plotters::coord::Shift>,
|
||
text: &str,
|
||
x: i32,
|
||
y: i32,
|
||
scale: i32,
|
||
color: &RGBColor,
|
||
) {
|
||
let mut cx = x;
|
||
for ch in text.chars() {
|
||
let glyph_ch = if BASIC_FONTS.get(ch).is_some() { ch } else { '?' };
|
||
if let Some(glyph) = BASIC_FONTS.get(glyph_ch) {
|
||
for (row, &bits) in glyph.iter().enumerate() {
|
||
for col in 0i32..8 {
|
||
if bits & (1 << col) != 0 {
|
||
for dy in 0..scale {
|
||
for dx in 0..scale {
|
||
area.draw_pixel(
|
||
(cx + col * scale + dx, y + row as i32 * scale + dy),
|
||
color,
|
||
)
|
||
.ok();
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
cx += 8 * scale;
|
||
}
|
||
}
|
||
|
||
/// Separable Gaussian blur on a flat f32 grid.
|
||
///
|
||
/// NaN pixels are treated as transparent: they contribute zero weight and the
|
||
/// result is renormalised over non-NaN neighbours only. Pixels with no non-NaN
|
||
/// neighbours remain NaN.
|
||
fn gaussian_blur(grid: &[f32], w: usize, h: usize, sigma: f32) -> Vec<f32> {
|
||
let radius = (3.0 * sigma).ceil() as usize;
|
||
let kernel: Vec<f32> = (0..=2 * radius)
|
||
.map(|i| {
|
||
let x = i as f32 - radius as f32;
|
||
(-x * x / (2.0 * sigma * sigma)).exp()
|
||
})
|
||
.collect();
|
||
|
||
// Horizontal pass
|
||
let temp: Vec<f32> = (0..h * w)
|
||
.into_par_iter()
|
||
.map(|idx| {
|
||
let row = idx / w;
|
||
let col = idx % w;
|
||
let mut val = 0.0f32;
|
||
let mut wsum = 0.0f32;
|
||
for (ki, &k) in kernel.iter().enumerate() {
|
||
let c = col as i64 + ki as i64 - radius as i64;
|
||
if c < 0 || c >= w as i64 { continue; }
|
||
let v = grid[row * w + c as usize];
|
||
if !v.is_nan() { val += k * v; wsum += k; }
|
||
}
|
||
if wsum > 0.0 { val / wsum } else { f32::NAN }
|
||
})
|
||
.collect();
|
||
|
||
// Vertical pass
|
||
(0..h * w)
|
||
.into_par_iter()
|
||
.map(|idx| {
|
||
let row = idx / w;
|
||
let col = idx % w;
|
||
let mut val = 0.0f32;
|
||
let mut wsum = 0.0f32;
|
||
for (ki, &k) in kernel.iter().enumerate() {
|
||
let r = row as i64 + ki as i64 - radius as i64;
|
||
if r < 0 || r >= h as i64 { continue; }
|
||
let v = temp[r as usize * w + col];
|
||
if !v.is_nan() { val += k * v; wsum += k; }
|
||
}
|
||
if wsum > 0.0 { val / wsum } else { f32::NAN }
|
||
})
|
||
.collect()
|
||
}
|
||
|
||
/// Render one forecast frame as a PNG.
|
||
///
|
||
/// - `lats`, `lons`, `cloud_cover` must all have the same length.
|
||
/// - `title` is displayed at the top of the image.
|
||
/// - Output is written to `output_path`.
|
||
pub fn render_frame(
|
||
output_path: &Path,
|
||
lats: &[f32],
|
||
lons: &[f32],
|
||
cloud_cover: &[f32],
|
||
title: &str,
|
||
) -> Result<()> {
|
||
let proj = OrthoProjection::new(CENTER_LAT, CENTER_LON, HALF_W_DEG, HALF_H_DEG);
|
||
|
||
let root = BitMapBackend::new(output_path, (IMG_WIDTH, IMG_HEIGHT)).into_drawing_area();
|
||
root.fill(&BACKGROUND).context("fill background")?;
|
||
|
||
const LEGEND_W: u32 = 130;
|
||
const TITLE_H: u32 = 40;
|
||
const MARGIN: u32 = 8;
|
||
|
||
let map_area = root.margin(TITLE_H, MARGIN, MARGIN, LEGEND_W + MARGIN);
|
||
|
||
let map_w = IMG_WIDTH - LEGEND_W - MARGIN * 2;
|
||
let map_h = IMG_HEIGHT - TITLE_H - MARGIN;
|
||
|
||
map_area.fill(&LAND_BG).context("fill map area")?;
|
||
|
||
// Build a Delaunay triangulation in continuous pixel-space coordinates.
|
||
// Include all grid points that project within a small margin outside the image
|
||
// bounds so that boundary pixels interpolate cleanly without edge artefacts.
|
||
const TRI_MARGIN: f64 = 10.0;
|
||
let mut tri = DelaunayTriangulation::<CloudVertex>::new();
|
||
for i in 0..lats.len() {
|
||
let cover = cloud_cover[i];
|
||
if cover.is_nan() {
|
||
continue;
|
||
}
|
||
let Some((px, py)) = proj.project(lats[i] as f64, lons[i] as f64) else { continue };
|
||
// Continuous pixel coordinates — no rounding.
|
||
let col_f = (px + HALF_W_DEG) / (2.0 * HALF_W_DEG) * map_w as f64;
|
||
let row_f = (-py + HALF_H_DEG) / (2.0 * HALF_H_DEG) * map_h as f64;
|
||
if col_f < -TRI_MARGIN
|
||
|| col_f > map_w as f64 + TRI_MARGIN
|
||
|| row_f < -TRI_MARGIN
|
||
|| row_f > map_h as f64 + TRI_MARGIN
|
||
{
|
||
continue;
|
||
}
|
||
let _ = tri.insert(CloudVertex { position: Point2::new(col_f, row_f), cloud: cover });
|
||
}
|
||
|
||
// For each output pixel, locate the containing triangle and barycentric-interpolate
|
||
// the cloud cover from its three vertices. The triangulation is immutable at this
|
||
// point, so the work is embarrassingly parallel via rayon.
|
||
let grid_w = map_w as usize;
|
||
let grid_h = map_h as usize;
|
||
let cover_grid: Vec<f32> = (0..grid_h * grid_w)
|
||
.into_par_iter()
|
||
.map(|i| {
|
||
let col = (i % grid_w) as f64 + 0.5;
|
||
let row = (i / grid_w) as f64 + 0.5;
|
||
let p = Point2::new(col, row);
|
||
use spade::PositionInTriangulation::*;
|
||
match tri.locate(p) {
|
||
OnVertex(v) => tri.vertex(v).data().cloud,
|
||
OnEdge(e) => {
|
||
// Linear interpolation along the edge.
|
||
let edge = tri.directed_edge(e);
|
||
let a = edge.from();
|
||
let b = edge.to();
|
||
let (ax, ay) = (a.position().x, a.position().y);
|
||
let (bx, by) = (b.position().x, b.position().y);
|
||
let ab2 = (bx - ax).powi(2) + (by - ay).powi(2);
|
||
let t = if ab2 == 0.0 {
|
||
0.0
|
||
} else {
|
||
(((p.x - ax) * (bx - ax) + (p.y - ay) * (by - ay)) / ab2)
|
||
.clamp(0.0, 1.0)
|
||
};
|
||
a.data().cloud * (1.0 - t) as f32 + b.data().cloud * t as f32
|
||
}
|
||
OnFace(f) => {
|
||
let face = tri.face(f);
|
||
let [v0, v1, v2] = face.vertices();
|
||
let (ax, ay) = (v0.position().x, v0.position().y);
|
||
let (bx, by) = (v1.position().x, v1.position().y);
|
||
let (cx, cy) = (v2.position().x, v2.position().y);
|
||
let denom = (by - cy) * (ax - cx) + (cx - bx) * (ay - cy);
|
||
if denom == 0.0 {
|
||
return v0.data().cloud;
|
||
}
|
||
let l0 = ((by - cy) * (p.x - cx) + (cx - bx) * (p.y - cy)) / denom;
|
||
let l1 = ((cy - ay) * (p.x - cx) + (ax - cx) * (p.y - cy)) / denom;
|
||
let l2 = 1.0 - l0 - l1;
|
||
v0.data().cloud * l0 as f32
|
||
+ v1.data().cloud * l1 as f32
|
||
+ v2.data().cloud * l2 as f32
|
||
}
|
||
OutsideOfConvexHull(_) | NoTriangulation => f32::NAN,
|
||
}
|
||
})
|
||
.collect();
|
||
|
||
const BLUR_SIGMA: f32 = 8.0;
|
||
let cover_grid = gaussian_blur(&cover_grid, grid_w, grid_h, BLUR_SIGMA);
|
||
|
||
// Paint the interpolated grid.
|
||
for row in 0..map_h as i32 {
|
||
for col in 0..map_w as i32 {
|
||
let cover = cover_grid[row as usize * grid_w + col as usize];
|
||
if let Some(color) = cloud_color(cover) {
|
||
map_area.draw_pixel((col, row), &color).ok();
|
||
}
|
||
// NaN or clear sky: keep the land background
|
||
}
|
||
}
|
||
|
||
// Draw city markers and labels
|
||
for &(label, lon, lat) in LOCATIONS {
|
||
let Some((px, py)) = proj.project(lat, lon) else { continue };
|
||
let Some((col, row)) = proj.to_pixel(px, py, map_w, map_h) else { continue };
|
||
|
||
// Small cross marker
|
||
for d in -3i32..=3i32 {
|
||
map_area.draw_pixel((col + d, row), &RED).ok();
|
||
map_area.draw_pixel((col, row + d), &RED).ok();
|
||
}
|
||
|
||
// Label offset slightly above and to the right of the marker
|
||
draw_pixel_text(&map_area, label, col + 5, row - 9, 1, &BLACK);
|
||
}
|
||
|
||
// Draw title (scale=2 → 16px tall)
|
||
draw_pixel_text(&root, title, MARGIN as i32, (TITLE_H / 2 - 8) as i32, 2, &BLACK);
|
||
|
||
// Draw legend
|
||
draw_legend(&root, IMG_WIDTH - LEGEND_W, TITLE_H)?;
|
||
|
||
root.present().context("write PNG")?;
|
||
Ok(())
|
||
}
|
||
|
||
fn draw_legend(
|
||
root: &DrawingArea<BitMapBackend, plotters::coord::Shift>,
|
||
x: u32,
|
||
y: u32,
|
||
) -> Result<()> {
|
||
draw_pixel_text(root, "Cloud cover", x as i32 + 4, y as i32 + 4, 1, &BLACK);
|
||
|
||
let box_h = 18i32;
|
||
let box_w = 24i32;
|
||
|
||
// Helper: draw one legend entry
|
||
let draw_entry = |label: &str, color: RGBColor, row_y: i32| {
|
||
root.draw(&Rectangle::new(
|
||
[(x as i32 + 4, row_y), (x as i32 + 4 + box_w, row_y + box_h)],
|
||
ShapeStyle { color: color.to_rgba(), filled: true, stroke_width: 0 },
|
||
))
|
||
.ok();
|
||
root.draw(&Rectangle::new(
|
||
[(x as i32 + 4, row_y), (x as i32 + 4 + box_w, row_y + box_h)],
|
||
ShapeStyle { color: BLACK.to_rgba(), filled: false, stroke_width: 1 },
|
||
))
|
||
.ok();
|
||
draw_pixel_text(
|
||
root,
|
||
label,
|
||
x as i32 + 4 + box_w + 4,
|
||
row_y + (box_h - 8) / 2,
|
||
1,
|
||
&BLACK,
|
||
);
|
||
};
|
||
|
||
let entries = [
|
||
(0.0f32, "< 1%"),
|
||
(25.0, "~25%"),
|
||
(50.0, "~50%"),
|
||
(75.0, "~75%"),
|
||
(100.0, "100%"),
|
||
];
|
||
|
||
for (i, &(cover, label)) in entries.iter().enumerate() {
|
||
let row_y = y as i32 + 22 + i as i32 * (box_h + 2);
|
||
let color = cloud_color(cover).unwrap_or(LAND_BG);
|
||
draw_entry(label, color, row_y);
|
||
}
|
||
|
||
Ok(())
|
||
}
|