schuwi/dwd_cloud_cover_rs
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Derive OSM viewport from zoom level instead of the other way around

Browse Source 
Replace the heuristic choose_zoom() back-computation with a fixed
OSM_ZOOM constant (zoom 10). The orthographic projection half-extents
are now derived so that one tile pixel maps exactly to one output pixel,
eliminating bilinear interpolation artefacts on the OSM base layer.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Schuwi
2026-03-07 11:04:26 +01:00
parent ff69c1d2db
commit 71d6588dc3
3 changed files with 29 additions and 35 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
src/main.rs
View File

@@ -134,16 +134,12 @@ fn main() -> Result<()> {
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 proj = render::make_projection();
let tile_cache = tiles::fetch_tiles(
&proj,
render::MAP_W,
render::MAP_H,
render::OSM_ZOOM,
&PathBuf::from(CACHE_DIR),
)?;
Some(tiles::rasterize_basemap(
@@ -157,16 +153,12 @@ fn main() -> Result<()> {
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 proj = render::make_projection();
let label_cache = tiles::fetch_label_tiles(
&proj,
render::MAP_W,
render::MAP_H,
render::OSM_ZOOM,
&PathBuf::from(CACHE_DIR),
)?;
Some(tiles::rasterize_labels(

31
src/render.rs
View File

@@ -15,13 +15,28 @@ 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.
/// Map image centre and OSM zoom level.
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);
/// OSM tile zoom level. The viewport half-extents are derived from this so that
/// one tile pixel maps to exactly one output pixel, eliminating interpolation
/// artefacts on the OSM base layer.
pub const OSM_ZOOM: u8 = 10;
/// Build the orthographic projection whose viewport is sized so that the
/// output resolution matches OSM tiles at `OSM_ZOOM` exactly.
///
/// At zoom `z`, Web-Mercator tiles cover `360 / 2^z` degrees of longitude per
/// tile, each 256 px wide. Near the centre latitude `φ`, 1° of longitude ≈
/// `cos φ` arc-degrees, so the arc-degree scale is:
/// `pixels_per_arc_deg = 256 · 2^z / (360 · cos φ)`
/// Half-extents follow directly from the map pixel dimensions.
pub fn make_projection() -> OrthoProjection {
let scale = f64::from(256u32 << OSM_ZOOM as u32) / (360.0 * CENTER_LAT.to_radians().cos());
let half_w = MAP_W as f64 / (2.0 * scale);
let half_h = MAP_H as f64 / (2.0 * scale);
OrthoProjection::new(CENTER_LAT, CENTER_LON, half_w, half_h)
}
const BACKGROUND: RGBColor = RGBColor(240, 248, 255); // pale sky-blue
const LAND_BG: RGBColor = RGBColor(142, 178, 35); // strong green for land
@@ -165,7 +180,7 @@ pub fn render_frame(
basemap: Option<&[[u8; 3]]>,
labels: Option<&[[u8; 4]]>,
) -> Result<()> {
let proj = OrthoProjection::new(CENTER_LAT, CENTER_LON, HALF_W_DEG, HALF_H_DEG);
let proj = make_projection();
let root = BitMapBackend::new(output_path, (IMG_WIDTH, IMG_HEIGHT)).into_drawing_area();
root.fill(&BACKGROUND).context("fill background")?;
@@ -188,8 +203,8 @@ pub fn render_frame(
}
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;
let col_f = (px + proj.half_width) / (2.0 * proj.half_width) * map_w as f64;
let row_f = (-py + proj.half_height) / (2.0 * proj.half_height) * map_h as f64;
if col_f < -TRI_MARGIN
|| col_f > map_w as f64 + TRI_MARGIN
|| row_f < -TRI_MARGIN

17
src/tiles.rs
View File

@@ -21,19 +21,6 @@ fn geo_to_tile_coords(lat_deg: f64, lon_deg: f64, zoom: u8) -> (f64, f64) {
(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)> {
@@ -171,10 +158,10 @@ pub fn fetch_tiles(
proj: &OrthoProjection,
map_w: u32,
map_h: u32,
zoom: u8,
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());
@@ -323,10 +310,10 @@ pub fn fetch_label_tiles(
proj: &OrthoProjection,
map_w: u32,
map_h: u32,
zoom: u8,
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());