Add OpenStreetMap base layer with Carto tile overlay
- Fetch Carto Voyager no-label base tiles and label-only tiles separately,
cached under cache/osm_tiles/{z}/ and cache/osm_tiles/labels/{z}/
- Rasterize both into the orthographic projection via exact inverse projection
(new unproject/pixel_to_geo methods on OrthoProjection)
- Render pipeline: basemap → cloud blend → label overlay (labels above clouds)
- Bilinear interpolation for base tiles; premultiplied-alpha bilinear for label
tiles to prevent dark-fringe artifacts at text edges
- Dynamic zoom selection (floor-based) from actual geographic bounding box
- Fix horizontal squish: derive HALF_W_DEG from pixel aspect ratio so
degrees-per-pixel is equal on both axes
- Add --no-basemap flag to skip tile fetching for offline/fast use
- Remove hardcoded city markers/labels when tile label overlay is present
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Schuwi
54
src/main.rs
54
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);
|
||||
|
||||
@@ -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.
|
||||
|
||||
101
src/render.rs
101
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)
|
||||
|
||||
367
src/tiles.rs
Normal file
367
src/tiles.rs
Normal file
@@ -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<RgbImage> {
|
||||
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<TileCache> {
|
||||
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<HashMap<_, _>> = 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<RgbaImage> {
|
||||
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<LabelTileCache> {
|
||||
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<HashMap<_, _>> = 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()
|
||||
}
|
||||
Reference in New Issue
Block a user