diff --git a/src/main.rs b/src/main.rs index c926bfc..8f845aa 100644 --- a/src/main.rs +++ b/src/main.rs @@ -2,6 +2,7 @@ mod download; mod grib; mod projection; mod render; +mod tiles; use anyhow::{bail, Context, Result}; use chrono::{DateTime, Datelike, Timelike, Utc}; @@ -29,6 +30,10 @@ struct Args { /// Number of forecast hours to download and render (0–48) #[arg(value_name = "HOURS")] prediction_hours: u32, + + /// Disable OpenStreetMap base layer (use flat green background) + #[arg(long)] + no_basemap: bool, } fn main() -> Result<()> { @@ -125,6 +130,53 @@ fn main() -> Result<()> { })?; cloud_data.sort_unstable_by_key(|(step, _)| *step); + // Fetch and rasterize OSM basemap (once, reused for all frames). + let basemap = if args.no_basemap { + None + } else { + let proj = projection::OrthoProjection::new( + render::CENTER_LAT, + render::CENTER_LON, + render::HALF_W_DEG, + render::HALF_H_DEG, + ); + let tile_cache = tiles::fetch_tiles( + &proj, + render::MAP_W, + render::MAP_H, + &PathBuf::from(CACHE_DIR), + )?; + Some(tiles::rasterize_basemap( + &proj, + render::MAP_W, + render::MAP_H, + &tile_cache, + )) + }; + + let labels = if args.no_basemap { + None + } else { + let proj = projection::OrthoProjection::new( + render::CENTER_LAT, + render::CENTER_LON, + render::HALF_W_DEG, + render::HALF_H_DEG, + ); + let label_cache = tiles::fetch_label_tiles( + &proj, + render::MAP_W, + render::MAP_H, + &PathBuf::from(CACHE_DIR), + )?; + Some(tiles::rasterize_labels( + &proj, + render::MAP_W, + render::MAP_H, + &label_cache, + )) + }; + eprintln!("Rendering {} frame(s)...", args.prediction_hours); let mut output_paths = Vec::new(); @@ -156,7 +208,7 @@ fn main() -> Result<()> { let out_name = format!("clct_{:06}.png", step + 1); let out_path = cache_dir.join(&out_name); - render::render_frame(&out_path, &lats, &lons, &cloud, &title) + render::render_frame(&out_path, &lats, &lons, &cloud, &title, basemap.as_deref(), labels.as_deref()) .with_context(|| format!("Rendering frame {}", step))?; eprintln!(" [{:02}/{:02}] {}", step + 1, args.prediction_hours, out_name); diff --git a/src/projection.rs b/src/projection.rs index fad9739..e20df33 100644 --- a/src/projection.rs +++ b/src/projection.rs @@ -50,6 +50,43 @@ impl OrthoProjection { Some((x.to_degrees(), y.to_degrees())) } + /// Inverse of `project`: convert projected (x, y) in degrees-of-arc back to + /// geographic (lat, lon) in degrees. + pub fn unproject(&self, x: f64, y: f64) -> (f64, f64) { + let x_rad = x.to_radians(); + let y_rad = y.to_radians(); + let rho = (x_rad * x_rad + y_rad * y_rad).sqrt(); + + if rho < 1e-14 { + // At the projection centre + return (self.phi0.to_degrees(), self.lam0.to_degrees()); + } + + let c = rho.asin(); + let sin_c = c.sin(); + let cos_c = c.cos(); + + let lat = (cos_c * self.phi0.sin() + y_rad * sin_c * self.phi0.cos() / rho).asin(); + let lon = self.lam0 + + (x_rad * sin_c).atan2(rho * self.phi0.cos() * cos_c - y_rad * self.phi0.sin() * sin_c); + + (lat.to_degrees(), lon.to_degrees()) + } + + /// Inverse of `to_pixel`: convert pixel (col, row) back to projected (x, y) + /// in degrees-of-arc. + pub fn from_pixel(&self, col: f64, row: f64, width: u32, height: u32) -> (f64, f64) { + let x = col / width as f64 * (2.0 * self.half_width) - self.half_width; + let y = self.half_height - row / height as f64 * (2.0 * self.half_height); + (x, y) + } + + /// Convenience: convert pixel (col, row) directly to geographic (lat, lon). + pub fn pixel_to_geo(&self, col: f64, row: f64, width: u32, height: u32) -> (f64, f64) { + let (x, y) = self.from_pixel(col, row, width, height); + self.unproject(x, y) + } + /// Convert projected (x, y) in degrees-of-arc to pixel (col, row). /// /// Returns `None` if the point falls outside the image boundary. diff --git a/src/render.rs b/src/render.rs index 50ed309..f97f935 100644 --- a/src/render.rs +++ b/src/render.rs @@ -6,14 +6,22 @@ use rayon::prelude::*; use spade::{DelaunayTriangulation, HasPosition, Point2, Triangulation}; use std::path::Path; -const IMG_WIDTH: u32 = 900; -const IMG_HEIGHT: u32 = 600; +pub const IMG_WIDTH: u32 = 900; +pub 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 +pub const LEGEND_W: u32 = 130; +pub const TITLE_H: u32 = 40; +pub const MARGIN: u32 = 8; +pub const MAP_W: u32 = IMG_WIDTH - LEGEND_W - MARGIN * 2; +pub const MAP_H: u32 = IMG_HEIGHT - TITLE_H - MARGIN; + +/// Map image centre and view window. +pub const CENTER_LAT: f64 = 52.56; +pub const CENTER_LON: f64 = 13.08; +pub const HALF_H_DEG: f64 = 0.4; // vertical half-extent in degrees of arc +// HALF_W_DEG is derived so that degrees-per-pixel is equal on both axes +// (MAP_W / MAP_H × HALF_H_DEG), preventing horizontal squish. +pub const HALF_W_DEG: f64 = HALF_H_DEG * (MAP_W as f64) / (MAP_H as f64); const BACKGROUND: RGBColor = RGBColor(240, 248, 255); // pale sky-blue const LAND_BG: RGBColor = RGBColor(142, 178, 35); // strong green for land @@ -154,22 +162,19 @@ pub fn render_frame( lons: &[f32], cloud_cover: &[f32], title: &str, + basemap: Option<&[[u8; 3]]>, + labels: Option<&[[u8; 4]]>, ) -> 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")?; + let map_w = MAP_W; + let map_h = MAP_H; + let n_pixels = map_w as usize * map_h as usize; // Build a Delaunay triangulation in continuous pixel-space coordinates. // Include all grid points that project within a small margin outside the image @@ -250,30 +255,58 @@ pub fn render_frame( 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 + // Build pixel buffer: start from basemap or flat LAND_BG. + let mut pixel_buf: Vec<[u8; 3]> = if let Some(bm) = basemap { + bm.to_vec() + } else { + vec![[LAND_BG.0, LAND_BG.1, LAND_BG.2]; n_pixels] + }; + + // Blend cloud cover toward white. + for (idx, &cover) in cover_grid.iter().enumerate() { + if cover.is_nan() || cover < 1.0 { + continue; + } + let t = (cover / 100.0).clamp(0.0, 1.0); + let [r, g, b] = pixel_buf[idx]; + let blend = |base: u8| -> u8 { (base as f32 + t * (255.0 - base as f32)).round() as u8 }; + pixel_buf[idx] = [blend(r), blend(g), blend(b)]; + } + + // Composite label overlay (map labels above clouds). + if let Some(lbl) = labels { + for (idx, &[lr, lg, lb, la]) in lbl.iter().enumerate() { + if la == 0 { continue; } + let a = la as f32 / 255.0; + let [r, g, b] = pixel_buf[idx]; + pixel_buf[idx] = [ + (lr as f32 * a + r as f32 * (1.0 - a)).round() as u8, + (lg as f32 * a + g as f32 * (1.0 - a)).round() as u8, + (lb as f32 * a + b as f32 * (1.0 - a)).round() as u8, + ]; } } - // 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(); + // Flush pixel buffer to drawing area. + for row in 0..map_h as i32 { + for col in 0..map_w as i32 { + let [r, g, b] = pixel_buf[row as usize * map_w as usize + col as usize]; + map_area.draw_pixel((col, row), &RGBColor(r, g, b)).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 city markers (cross only; labels come from the tile label overlay) + if labels.is_none() { + 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 }; + + for d in -3i32..=3i32 { + map_area.draw_pixel((col + d, row), &RED).ok(); + map_area.draw_pixel((col, row + d), &RED).ok(); + } + draw_pixel_text(&map_area, label, col + 5, row - 9, 1, &BLACK); + } } // Draw title (scale=2 → 16px tall) diff --git a/src/tiles.rs b/src/tiles.rs new file mode 100644 index 0000000..65cc84a --- /dev/null +++ b/src/tiles.rs @@ -0,0 +1,367 @@ +use crate::projection::OrthoProjection; +use anyhow::{bail, Context, Result}; +use image::{RgbaImage, RgbImage}; +use rayon::prelude::*; +use std::collections::HashMap; +use std::path::Path; +use std::sync::Mutex; + +const TILE_SIZE: u32 = 256; +const TILE_URL_BASE: &str = "https://cartodb-basemaps-a.global.ssl.fastly.net/rastertiles/voyager_nolabels"; +const TILE_URL_LABELS: &str = "https://cartodb-basemaps-a.global.ssl.fastly.net/rastertiles/voyager_only_labels"; +const USER_AGENT: &str = "cloud_cover_rs/0.1 (DWD cloud cover visualization)"; + +/// Convert geographic (lat, lon) in degrees to fractional tile coordinates at +/// the given zoom level (slippy map / Web Mercator convention). +fn geo_to_tile_coords(lat_deg: f64, lon_deg: f64, zoom: u8) -> (f64, f64) { + let n = f64::from(1u32 << zoom); + let x = (lon_deg + 180.0) / 360.0 * n; + let lat_rad = lat_deg.to_radians(); + let y = (1.0 - (lat_rad.tan() + 1.0 / lat_rad.cos()).ln() / std::f64::consts::PI) / 2.0 * n; + (x, y) +} + +/// Choose a zoom level so that tile pixels approximately match the output +/// resolution. `lon_range` is the geographic longitude span in degrees, +/// `map_w` is the output width in pixels. +fn choose_zoom(lon_range: f64, map_w: u32) -> u8 { + // At zoom z, the world is 256 * 2^z pixels wide covering 360° of longitude. + // We want: 256 * 2^z / 360 >= map_w / lon_range (tiles at least as detailed as output) + // => 2^z >= map_w * 360 / (lon_range * 256) + // Use floor so we pick the level just below 1:1, then bilinear upscaling gives clean result. + let needed = (map_w as f64 * 360.0) / (lon_range * TILE_SIZE as f64); + let z = needed.log2().floor() as u8; + z.clamp(1, 18) +} + +/// Determine the set of tile (x, y) indices needed to cover the given +/// geographic bounding box at the given zoom level. +fn required_tiles(min_lon: f64, max_lon: f64, min_lat: f64, max_lat: f64, zoom: u8) -> Vec<(u32, u32)> { + // Note: higher latitude → lower y-tile (north = top) + let (x_min_f, y_min_f) = geo_to_tile_coords(max_lat, min_lon, zoom); + let (x_max_f, y_max_f) = geo_to_tile_coords(min_lat, max_lon, zoom); + + let x_min = x_min_f.floor() as u32; + let x_max = x_max_f.floor() as u32; + let y_min = y_min_f.floor() as u32; + let y_max = y_max_f.floor() as u32; + + let mut tiles = Vec::new(); + for x in x_min..=x_max { + for y in y_min..=y_max { + tiles.push((x, y)); + } + } + tiles +} + +/// Compute the geographic bounding box for the map area defined by `proj`. +fn map_bbox(proj: &OrthoProjection, map_w: u32, map_h: u32) -> (f64, f64, f64, f64) { + let w = map_w as f64; + let h = map_h as f64; + let mut min_lat = f64::MAX; + let mut max_lat = f64::MIN; + let mut min_lon = f64::MAX; + let mut max_lon = f64::MIN; + for (col, row) in [ + (0.0, 0.0), (w, 0.0), (0.0, h), (w, h), + (w / 2.0, 0.0), (w / 2.0, h), + (0.0, h / 2.0), (w, h / 2.0), + ] { + let (lat, lon) = proj.pixel_to_geo(col, row, map_w, map_h); + min_lat = min_lat.min(lat); + max_lat = max_lat.max(lat); + min_lon = min_lon.min(lon); + max_lon = max_lon.max(lon); + } + (min_lon, max_lon, min_lat, max_lat) +} + +// ─── RGB tile cache (base layer, no labels) ────────────────────────────────── + +/// A cache of decoded RGB tiles (no labels), keyed by (x_tile, y_tile). +pub struct TileCache { + zoom: u8, + tiles: HashMap<(u32, u32), RgbImage>, +} + +impl TileCache { + /// Sample the tile cache at a geographic coordinate using bilinear interpolation. + pub fn sample(&self, lat: f64, lon: f64) -> Option<[u8; 3]> { + let (tx, ty) = geo_to_tile_coords(lat, lon, self.zoom); + + // Absolute sub-pixel position in the tile grid (in tile pixels) + let apx = tx * TILE_SIZE as f64; + let apy = ty * TILE_SIZE as f64; + let px0 = apx.floor() as u32; + let py0 = apy.floor() as u32; + let px1 = px0 + 1; + let py1 = py0 + 1; + let fx = (apx - px0 as f64) as f32; + let fy = (apy - py0 as f64) as f32; + + let get = |apx: u32, apy: u32| -> Option<[u8; 3]> { + let tx = apx / TILE_SIZE; + let ty = apy / TILE_SIZE; + let lx = apx % TILE_SIZE; + let ly = apy % TILE_SIZE; + let tile = self.tiles.get(&(tx, ty))?; + Some(tile.get_pixel(lx, ly).0) + }; + + let p00 = get(px0, py0)?; + let p10 = get(px1, py0).unwrap_or(p00); + let p01 = get(px0, py1).unwrap_or(p00); + let p11 = get(px1, py1).unwrap_or(p00); + + let blend_ch = |c00: u8, c10: u8, c01: u8, c11: u8| -> u8 { + let top = c00 as f32 * (1.0 - fx) + c10 as f32 * fx; + let bot = c01 as f32 * (1.0 - fx) + c11 as f32 * fx; + (top * (1.0 - fy) + bot * fy).round() as u8 + }; + + Some([ + blend_ch(p00[0], p10[0], p01[0], p11[0]), + blend_ch(p00[1], p10[1], p01[1], p11[1]), + blend_ch(p00[2], p10[2], p01[2], p11[2]), + ]) + } +} + +fn fetch_rgb_tile(x: u32, y: u32, zoom: u8, url_base: &str, cache_dir: &Path) -> Result { + let cache_path = cache_dir + .join(zoom.to_string()) + .join(x.to_string()) + .join(format!("{y}.png")); + + if cache_path.exists() { + return Ok(image::open(&cache_path) + .with_context(|| format!("Failed to decode cached tile {}", cache_path.display()))? + .into_rgb8()); + } + + let url = format!("{url_base}/{zoom}/{x}/{y}.png"); + eprintln!(" Fetching tile {url}"); + + let client = reqwest::blocking::Client::builder() + .user_agent(USER_AGENT) + .build() + .context("Failed to build HTTP client")?; + let response = client.get(&url).send() + .with_context(|| format!("HTTP request failed for {url}"))?; + let status = response.status(); + if !status.is_success() { + bail!("HTTP {status} for {url}"); + } + let bytes = response.bytes() + .with_context(|| format!("Failed to read response body from {url}"))?; + + if let Some(parent) = cache_path.parent() { + std::fs::create_dir_all(parent) + .with_context(|| format!("Failed to create tile cache dir {}", parent.display()))?; + } + std::fs::write(&cache_path, &bytes) + .with_context(|| format!("Failed to write tile cache {}", cache_path.display()))?; + + Ok(image::load_from_memory(&bytes).context("Failed to decode tile PNG")?.into_rgb8()) +} + +/// Fetch all no-label base tiles needed to cover the map area. +pub fn fetch_tiles( + proj: &OrthoProjection, + map_w: u32, + map_h: u32, + cache_dir: &Path, +) -> Result { + let (min_lon, max_lon, min_lat, max_lat) = map_bbox(proj, map_w, map_h); + let zoom = choose_zoom(max_lon - min_lon, map_w); + let tile_coords = required_tiles(min_lon, max_lon, min_lat, max_lat, zoom); + eprintln!("Fetching {} base tiles (zoom {zoom})...", tile_coords.len()); + + let tile_cache_dir = cache_dir.join("osm_tiles"); + let tiles: Mutex> = Mutex::new(HashMap::new()); + + let pool = rayon::ThreadPoolBuilder::new().num_threads(2).build() + .context("Failed to build tile fetch thread pool")?; + pool.install(|| { + tile_coords.par_iter().try_for_each(|&(x, y)| -> Result<()> { + let img = fetch_rgb_tile(x, y, zoom, TILE_URL_BASE, &tile_cache_dir)?; + tiles.lock().unwrap().insert((x, y), img); + Ok(()) + }) + })?; + + Ok(TileCache { zoom, tiles: tiles.into_inner().unwrap() }) +} + +/// Rasterize the no-label basemap for the entire map area. +pub fn rasterize_basemap( + proj: &OrthoProjection, + map_w: u32, + map_h: u32, + tile_cache: &TileCache, +) -> Vec<[u8; 3]> { + let size = map_w as usize * map_h as usize; + (0..size) + .into_par_iter() + .map(|i| { + let col = (i % map_w as usize) as f64 + 0.5; + let row = (i / map_w as usize) as f64 + 0.5; + let (lat, lon) = proj.pixel_to_geo(col, row, map_w, map_h); + tile_cache.sample(lat, lon).unwrap_or([200, 200, 200]) + }) + .collect() +} + +// ─── RGBA label tile cache (labels only, transparent background) ────────────── + +/// A cache of decoded RGBA label-only tiles, keyed by (x_tile, y_tile). +pub struct LabelTileCache { + zoom: u8, + tiles: HashMap<(u32, u32), RgbaImage>, +} + +impl LabelTileCache { + /// Sample the label cache at a geographic coordinate using bilinear interpolation + /// in premultiplied-alpha space, then convert back to straight alpha. + /// + /// Premultiplied interpolation prevents dark fringe artifacts: transparent pixels + /// with undefined RGB values (commonly [0,0,0,0]) would otherwise bleed dark + /// colours into semi-transparent text edges in straight-alpha blending. + pub fn sample(&self, lat: f64, lon: f64) -> Option<[u8; 4]> { + let (tx, ty) = geo_to_tile_coords(lat, lon, self.zoom); + + let apx = tx * TILE_SIZE as f64; + let apy = ty * TILE_SIZE as f64; + let px0 = apx.floor() as u32; + let py0 = apy.floor() as u32; + let px1 = px0 + 1; + let py1 = py0 + 1; + let fx = (apx - px0 as f64) as f32; + let fy = (apy - py0 as f64) as f32; + + let get = |apx: u32, apy: u32| -> Option<[f32; 4]> { + let tx = apx / TILE_SIZE; + let ty = apy / TILE_SIZE; + let lx = apx % TILE_SIZE; + let ly = apy % TILE_SIZE; + let tile = self.tiles.get(&(tx, ty))?; + let [r, g, b, a] = tile.get_pixel(lx, ly).0; + // Convert to premultiplied float + let af = a as f32 / 255.0; + Some([r as f32 * af, g as f32 * af, b as f32 * af, a as f32]) + }; + + let p00 = get(px0, py0)?; + let p10 = get(px1, py0).unwrap_or(p00); + let p01 = get(px0, py1).unwrap_or(p00); + let p11 = get(px1, py1).unwrap_or(p00); + + let blend_ch = |c00: f32, c10: f32, c01: f32, c11: f32| -> f32 { + let top = c00 * (1.0 - fx) + c10 * fx; + let bot = c01 * (1.0 - fx) + c11 * fx; + top * (1.0 - fy) + bot * fy + }; + + let pr = blend_ch(p00[0], p10[0], p01[0], p11[0]); + let pg = blend_ch(p00[1], p10[1], p01[1], p11[1]); + let pb = blend_ch(p00[2], p10[2], p01[2], p11[2]); + let pa = blend_ch(p00[3], p10[3], p01[3], p11[3]); + + // Convert back from premultiplied to straight alpha + let (r, g, b) = if pa > 0.0 { + let inv = 255.0 / pa; + ((pr * inv).round() as u8, (pg * inv).round() as u8, (pb * inv).round() as u8) + } else { + (0, 0, 0) + }; + + Some([r, g, b, pa.round() as u8]) + } +} + +fn fetch_rgba_tile(x: u32, y: u32, zoom: u8, url_base: &str, cache_dir: &Path) -> Result { + let cache_path = cache_dir + .join(zoom.to_string()) + .join(x.to_string()) + .join(format!("{y}.png")); + + if cache_path.exists() { + return Ok(image::open(&cache_path) + .with_context(|| format!("Failed to decode cached label tile {}", cache_path.display()))? + .into_rgba8()); + } + + let url = format!("{url_base}/{zoom}/{x}/{y}.png"); + eprintln!(" Fetching label tile {url}"); + + let client = reqwest::blocking::Client::builder() + .user_agent(USER_AGENT) + .build() + .context("Failed to build HTTP client")?; + let response = client.get(&url).send() + .with_context(|| format!("HTTP request failed for {url}"))?; + let status = response.status(); + if !status.is_success() { + bail!("HTTP {status} for {url}"); + } + let bytes = response.bytes() + .with_context(|| format!("Failed to read response body from {url}"))?; + + if let Some(parent) = cache_path.parent() { + std::fs::create_dir_all(parent) + .with_context(|| format!("Failed to create label tile cache dir {}", parent.display()))?; + } + std::fs::write(&cache_path, &bytes) + .with_context(|| format!("Failed to write label tile cache {}", cache_path.display()))?; + + Ok(image::load_from_memory(&bytes).context("Failed to decode label tile PNG")?.into_rgba8()) +} + +/// Fetch all label-only tiles needed to cover the map area. +pub fn fetch_label_tiles( + proj: &OrthoProjection, + map_w: u32, + map_h: u32, + cache_dir: &Path, +) -> Result { + let (min_lon, max_lon, min_lat, max_lat) = map_bbox(proj, map_w, map_h); + let zoom = choose_zoom(max_lon - min_lon, map_w); + let tile_coords = required_tiles(min_lon, max_lon, min_lat, max_lat, zoom); + eprintln!("Fetching {} label tiles (zoom {zoom})...", tile_coords.len()); + + let tile_cache_dir = cache_dir.join("osm_tiles").join("labels"); + let tiles: Mutex> = Mutex::new(HashMap::new()); + + let pool = rayon::ThreadPoolBuilder::new().num_threads(2).build() + .context("Failed to build tile fetch thread pool")?; + pool.install(|| { + tile_coords.par_iter().try_for_each(|&(x, y)| -> Result<()> { + let img = fetch_rgba_tile(x, y, zoom, TILE_URL_LABELS, &tile_cache_dir)?; + tiles.lock().unwrap().insert((x, y), img); + Ok(()) + }) + })?; + + Ok(LabelTileCache { zoom, tiles: tiles.into_inner().unwrap() }) +} + +/// Rasterize the label overlay for the entire map area. +/// Returns one RGBA per pixel (alpha=0 where there is no label). +pub fn rasterize_labels( + proj: &OrthoProjection, + map_w: u32, + map_h: u32, + tile_cache: &LabelTileCache, +) -> Vec<[u8; 4]> { + let size = map_w as usize * map_h as usize; + (0..size) + .into_par_iter() + .map(|i| { + let col = (i % map_w as usize) as f64 + 0.5; + let row = (i / map_w as usize) as f64 + 0.5; + let (lat, lon) = proj.pixel_to_geo(col, row, map_w, map_h); + tile_cache.sample(lat, lon).unwrap_or([0, 0, 0, 0]) + }) + .collect() +}