Compare commits
11 Commits
6c9a20bf59
...
56f779990d
| Author | SHA1 | Date | |
|---|---|---|---|
|
56f779990d | ||
|
|
28b7e90af8 | ||
|
|
c7c08ebead | ||
|
|
71d6588dc3 | ||
|
|
ff69c1d2db | ||
|
|
ac9f9da3c4 | ||
|
|
764f4f6378 | ||
|
|
0c73e04959 | ||
|
|
ee9b039d56 | ||
|
|
4ee33a882d | ||
|
|
a9eed7f378 |
205
Cargo.lock
generated
205
Cargo.lock
generated
@@ -8,6 +8,12 @@ version = "2.0.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa"
|
||||
|
||||
[[package]]
|
||||
name = "allocator-api2"
|
||||
version = "0.2.21"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "683d7910e743518b0e34f1186f92494becacb047c7b6bf616c96772180fef923"
|
||||
|
||||
[[package]]
|
||||
name = "android_system_properties"
|
||||
version = "0.1.5"
|
||||
@@ -121,6 +127,12 @@ version = "1.5.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b"
|
||||
|
||||
[[package]]
|
||||
name = "byteorder-lite"
|
||||
version = "0.1.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "8f1fe948ff07f4bd06c30984e69f5b4899c516a3ef74f34df92a2df2ab535495"
|
||||
|
||||
[[package]]
|
||||
name = "bytes"
|
||||
version = "1.11.1"
|
||||
@@ -180,6 +192,16 @@ dependencies = [
|
||||
"windows-link",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "chrono-tz"
|
||||
version = "0.10.4"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "a6139a8597ed92cf816dfb33f5dd6cf0bb93a6adc938f11039f371bc5bcd26c3"
|
||||
dependencies = [
|
||||
"chrono",
|
||||
"phf",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "clap"
|
||||
version = "4.5.60"
|
||||
@@ -227,10 +249,15 @@ dependencies = [
|
||||
"anyhow",
|
||||
"bzip2",
|
||||
"chrono",
|
||||
"chrono-tz",
|
||||
"clap",
|
||||
"font8x8",
|
||||
"gif",
|
||||
"image 0.25.9",
|
||||
"plotters",
|
||||
"rayon",
|
||||
"reqwest",
|
||||
"spade",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
@@ -260,6 +287,31 @@ dependencies = [
|
||||
"cfg-if",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-deque"
|
||||
version = "0.8.6"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51"
|
||||
dependencies = [
|
||||
"crossbeam-epoch",
|
||||
"crossbeam-utils",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-epoch"
|
||||
version = "0.9.18"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e"
|
||||
dependencies = [
|
||||
"crossbeam-utils",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "crossbeam-utils"
|
||||
version = "0.8.21"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28"
|
||||
|
||||
[[package]]
|
||||
name = "displaydoc"
|
||||
version = "0.2.5"
|
||||
@@ -271,6 +323,18 @@ dependencies = [
|
||||
"syn",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "either"
|
||||
version = "1.15.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719"
|
||||
|
||||
[[package]]
|
||||
name = "equivalent"
|
||||
version = "1.0.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "877a4ace8713b0bcf2a4e7eec82529c029f1d0619886d18145fea96c3ffe5c0f"
|
||||
|
||||
[[package]]
|
||||
name = "fdeflate"
|
||||
version = "0.3.7"
|
||||
@@ -296,6 +360,12 @@ dependencies = [
|
||||
"miniz_oxide",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "foldhash"
|
||||
version = "0.1.5"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "d9c4f5dac5e15c24eb999c26181a6ca40b39fe946cbe4c263c7209467bc83af2"
|
||||
|
||||
[[package]]
|
||||
name = "font8x8"
|
||||
version = "0.2.7"
|
||||
@@ -387,6 +457,27 @@ dependencies = [
|
||||
"wasm-bindgen",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "gif"
|
||||
version = "0.14.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "f5df2ba84018d80c213569363bdcd0c64e6933c67fe4c1d60ecf822971a3c35e"
|
||||
dependencies = [
|
||||
"color_quant",
|
||||
"weezl",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "hashbrown"
|
||||
version = "0.15.5"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "9229cfe53dfd69f0609a49f65461bd93001ea1ef889cd5529dd176593f5338a1"
|
||||
dependencies = [
|
||||
"allocator-api2",
|
||||
"equivalent",
|
||||
"foldhash",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "heck"
|
||||
version = "0.5.0"
|
||||
@@ -630,7 +721,20 @@ dependencies = [
|
||||
"color_quant",
|
||||
"jpeg-decoder",
|
||||
"num-traits",
|
||||
"png",
|
||||
"png 0.17.16",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "image"
|
||||
version = "0.25.9"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "e6506c6c10786659413faa717ceebcb8f70731c0a60cbae39795fdf114519c1a"
|
||||
dependencies = [
|
||||
"bytemuck",
|
||||
"byteorder-lite",
|
||||
"moxcms",
|
||||
"num-traits",
|
||||
"png 0.18.1",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
@@ -728,6 +832,16 @@ dependencies = [
|
||||
"windows-sys 0.61.2",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "moxcms"
|
||||
version = "0.7.11"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "ac9557c559cd6fc9867e122e20d2cbefc9ca29d80d027a8e39310920ed2f0a97"
|
||||
dependencies = [
|
||||
"num-traits",
|
||||
"pxfm",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "num-traits"
|
||||
version = "0.2.19"
|
||||
@@ -755,6 +869,24 @@ version = "2.3.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "9b4f627cb1b25917193a259e49bdad08f671f8d9708acfd5fe0a8c1455d87220"
|
||||
|
||||
[[package]]
|
||||
name = "phf"
|
||||
version = "0.12.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "913273894cec178f401a31ec4b656318d95473527be05c0752cc41cdc32be8b7"
|
||||
dependencies = [
|
||||
"phf_shared",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "phf_shared"
|
||||
version = "0.12.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "06005508882fb681fd97892ecff4b7fd0fee13ef1aa569f8695dae7ab9099981"
|
||||
dependencies = [
|
||||
"siphasher",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "pin-project-lite"
|
||||
version = "0.2.17"
|
||||
@@ -798,7 +930,7 @@ version = "0.3.7"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "72ce181e3f6bf82d6c1dc569103ca7b1bd964c60ba03d7e6cdfbb3e3eb7f7405"
|
||||
dependencies = [
|
||||
"image",
|
||||
"image 0.24.9",
|
||||
"plotters-backend",
|
||||
]
|
||||
|
||||
@@ -815,6 +947,19 @@ dependencies = [
|
||||
"miniz_oxide",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "png"
|
||||
version = "0.18.1"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "60769b8b31b2a9f263dae2776c37b1b28ae246943cf719eb6946a1db05128a61"
|
||||
dependencies = [
|
||||
"bitflags 2.11.0",
|
||||
"crc32fast",
|
||||
"fdeflate",
|
||||
"flate2",
|
||||
"miniz_oxide",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "potential_utf"
|
||||
version = "0.1.4"
|
||||
@@ -842,6 +987,12 @@ dependencies = [
|
||||
"unicode-ident",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "pxfm"
|
||||
version = "0.1.28"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "b5a041e753da8b807c9255f28de81879c78c876392ff2469cde94799b2896b9d"
|
||||
|
||||
[[package]]
|
||||
name = "quinn"
|
||||
version = "0.11.9"
|
||||
@@ -941,6 +1092,26 @@ dependencies = [
|
||||
"getrandom 0.3.4",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rayon"
|
||||
version = "1.11.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "368f01d005bf8fd9b1206fb6fa653e6c4a81ceb1466406b81792d87c5677a58f"
|
||||
dependencies = [
|
||||
"either",
|
||||
"rayon-core",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "rayon-core"
|
||||
version = "1.13.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91"
|
||||
dependencies = [
|
||||
"crossbeam-deque",
|
||||
"crossbeam-utils",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "reqwest"
|
||||
version = "0.12.28"
|
||||
@@ -995,6 +1166,12 @@ dependencies = [
|
||||
"windows-sys 0.52.0",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "robust"
|
||||
version = "1.2.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "4e27ee8bb91ca0adcf0ecb116293afa12d393f9c2b9b9cd54d33e8078fe19839"
|
||||
|
||||
[[package]]
|
||||
name = "rustc-hash"
|
||||
version = "2.1.1"
|
||||
@@ -1115,6 +1292,12 @@ version = "0.3.8"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "e320a6c5ad31d271ad523dcf3ad13e2767ad8b1cb8f047f75a8aeaf8da139da2"
|
||||
|
||||
[[package]]
|
||||
name = "siphasher"
|
||||
version = "1.0.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "b2aa850e253778c88a04c3d7323b043aeda9d3e30d5971937c1855769763678e"
|
||||
|
||||
[[package]]
|
||||
name = "slab"
|
||||
version = "0.4.12"
|
||||
@@ -1137,6 +1320,18 @@ dependencies = [
|
||||
"windows-sys 0.60.2",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "spade"
|
||||
version = "2.15.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "fb313e1c8afee5b5647e00ee0fe6855e3d529eb863a0fdae1d60006c4d1e9990"
|
||||
dependencies = [
|
||||
"hashbrown",
|
||||
"num-traits",
|
||||
"robust",
|
||||
"smallvec",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "stable_deref_trait"
|
||||
version = "1.2.1"
|
||||
@@ -1473,6 +1668,12 @@ dependencies = [
|
||||
"rustls-pki-types",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "weezl"
|
||||
version = "0.1.12"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "a28ac98ddc8b9274cb41bb4d9d4d5c425b6020c50c46f25559911905610b4a88"
|
||||
|
||||
[[package]]
|
||||
name = "windows-core"
|
||||
version = "0.62.2"
|
||||
|
||||
@@ -15,3 +15,8 @@ font8x8 = "0.2"
|
||||
clap = { version = "4", features = ["derive"] }
|
||||
chrono = { version = "0.4", default-features = false, features = ["std", "clock"] }
|
||||
anyhow = "1"
|
||||
spade = "2"
|
||||
rayon = "1"
|
||||
image = { version = "0.25", default-features = false, features = ["png"] }
|
||||
gif = "0.14"
|
||||
chrono-tz = "0.10.4"
|
||||
|
||||
156
README.md
156
README.md
@@ -1,97 +1,141 @@
|
||||
# cloud_cover
|
||||
|
||||
Downloads total cloud cover (CLCT) predictions from the [DWD open-data
|
||||
server](https://opendata.dwd.de/weather/nwp/icon-d2/grib/) and renders
|
||||
one PNG map per forecast hour, centred on the Brandenburg / Berlin area.
|
||||
Cloud cover forecast maps for astrophotography planning, powered by DWD open data.
|
||||
|
||||
This is a Rust port of `../cloud_cover_prediction.py`. The output format
|
||||
differs: instead of using PyNGL's contour renderer the Rust version does a
|
||||
direct point-by-point orthographic projection of the ICON-D2 grid onto a
|
||||
bitmap, which avoids every system-library dependency (no CDO, no PyNIO, no
|
||||
PyNGL, no libssl, no libbz2 — the binary is fully self-contained).
|
||||
|
||||
|
||||
## What is this?
|
||||
|
||||
The German Weather Service (DWD) publishes free, high-resolution numerical weather
|
||||
predictions through its [open-data server](https://opendata.dwd.de/). The
|
||||
**ICON-D2** model covers Germany and surrounding areas with ~2 km grid spacing and
|
||||
produces forecasts up to 48 hours ahead, updated every 3 hours.
|
||||
|
||||
This tool downloads the **total cloud cover (CLCT)** field from ICON-D2 and renders
|
||||
it as a series of map images — one per forecast hour — overlaid on an OpenStreetMap
|
||||
base layer. The result is a quick visual answer to "will the sky be clear tonight?"
|
||||
|
||||
## Quick start
|
||||
|
||||
```sh
|
||||
cargo build --release
|
||||
./target/release/cloud_cover "2026-03-07T09:00:00Z" 24
|
||||
```
|
||||
|
||||
Output appears in `cache/2026-03-07_09UTC/` — one PNG per forecast hour plus an
|
||||
animated GIF.
|
||||
|
||||
## Usage
|
||||
|
||||
```
|
||||
cloud_cover <ISO_TIMESTAMP> <HOURS>
|
||||
cloud_cover [OPTIONS] <ISO_TIMESTAMP> <HOURS>
|
||||
```
|
||||
|
||||
| Argument | Description |
|
||||
|---|---|
|
||||
| `ISO_TIMESTAMP` | Model run time in ISO 8601, e.g. `2026-03-03T18:00:00Z` |
|
||||
| `ISO_TIMESTAMP` | Model run time in ISO 8601, e.g. `2026-03-07T09:00:00Z` |
|
||||
| `HOURS` | Number of forecast steps to render (0–48) |
|
||||
|
||||
The timestamp must correspond to an ICON-D2 model run, which is issued
|
||||
every 3 hours starting at 00 UTC (00, 03, 06, 09, 12, 15, 18, 21).
|
||||
The timestamp must match an ICON-D2 model run: every 3 hours starting at 00 UTC
|
||||
(00, 03, 06, 09, 12, 15, 18, 21).
|
||||
|
||||
### Example
|
||||
### Options
|
||||
|
||||
| Flag | Default | Description |
|
||||
|---|---|---|
|
||||
| `--center-lat` | 52.56 | Map centre latitude (°N) |
|
||||
| `--center-lon` | 13.08 | Map centre longitude (°E) |
|
||||
| `--zoom` | 10 | OSM tile zoom level (higher = more detail, smaller area) |
|
||||
| `--no-basemap` | — | Skip OSM tiles; use a plain green background |
|
||||
|
||||
### Examples
|
||||
|
||||
```sh
|
||||
cargo run --release -- "2026-03-03T18:00:00Z" 12
|
||||
# Next 12 hours from the 18 UTC run, default viewport (Falkensee / Berlin)
|
||||
./target/release/cloud_cover "2026-03-06T18:00:00Z" 12
|
||||
|
||||
# Wider area centred on central Germany at zoom 8
|
||||
./target/release/cloud_cover "2026-03-07T06:00:00Z" 24 \
|
||||
--center-lat 51.0 --center-lon 10.5 --zoom 8
|
||||
|
||||
# Quick run without map tiles
|
||||
./target/release/cloud_cover "2026-03-07T09:00:00Z" 6 --no-basemap
|
||||
```
|
||||
|
||||
Output PNGs are written to `cache/<date>_<hour>UTC/clct_000001.png` …
|
||||
`clct_NNNNNN.png`. Downloaded GRIB2 files are cached in the same directory
|
||||
so that re-running the same timestamp skips the network requests.
|
||||
## Output
|
||||
|
||||
Each run produces files in `cache/<date>_<hour>UTC/`:
|
||||
|
||||
- **`clct_000001.png` … `clct_NNNNNN.png`** — one 900 × 600 px map per forecast
|
||||
hour, showing cloud cover blended over the OSM base layer with a colour legend
|
||||
- **`clct_animation.gif`** — animated loop of all frames (generated when there are
|
||||
multiple steps)
|
||||
- **`clat.grib2`, `clon.grib2`, `clct_NNN.grib2`** — cached raw forecast data;
|
||||
re-running the same timestamp skips the download
|
||||
|
||||
Image titles show the forecast time converted to CET/CEST (Europe/Berlin). The first
|
||||
frame is labelled "Conditions at …", subsequent frames "Prediction at … (+NNh)".
|
||||
|
||||
Cloud cover is rendered as a continuous blend from the base map (clear sky) toward a
|
||||
light blue-white tone (overcast), matching the style familiar from weather apps.
|
||||
|
||||
## How it works
|
||||
|
||||
1. **Download** — `reqwest` fetches the bzip2-compressed GRIB2 files from
|
||||
the DWD server. Two coordinate files (`clat`, `clon`) give the latitude
|
||||
and longitude of every point on the ICON-D2 icosahedral grid (~542 000
|
||||
points for Germany). One `clct` file per forecast step contains the
|
||||
total cloud cover percentage at each point.
|
||||
1. **Download** — Fetches bzip2-compressed GRIB2 files from the DWD open-data
|
||||
server: two coordinate grids (`clat`, `clon`) describing the ~542,000 points of
|
||||
the ICON-D2 icosahedral grid over Germany, plus one `clct` file per forecast step.
|
||||
Downloads run in parallel and are cached to disk.
|
||||
|
||||
2. **Decompress** — `bzip2` decompresses the files in-memory.
|
||||
2. **Parse** — A built-in minimal GRIB2 decoder extracts the data arrays. Only
|
||||
simple packing (data representation template 0) is implemented, which is the
|
||||
format DWD uses for these fields. Grid points flagged absent by the GRIB2 bitmap
|
||||
are set to NaN.
|
||||
|
||||
3. **Parse** — A minimal built-in GRIB2 decoder extracts the data values.
|
||||
Only simple packing (data representation template 0) is supported, which
|
||||
is what DWD uses. Grid points that are flagged as absent by the section-6
|
||||
bitmap are set to `NaN` and skipped during rendering.
|
||||
3. **Interpolate** — The irregular grid points are projected into pixel space via an
|
||||
orthographic projection, then connected into a Delaunay triangulation. Each output
|
||||
pixel is interpolated via barycentric coordinates within its enclosing triangle,
|
||||
followed by a NaN-aware Gaussian blur to smooth triangle edges.
|
||||
|
||||
4. **Project & render** — Each grid point is projected onto the image plane
|
||||
using an orthographic projection centred on Falkensee (52.56 °N,
|
||||
13.08 °E). Points are painted as 2 × 2 pixel squares, colour-coded by
|
||||
cloud cover percentage. City markers and labels are drawn on top, and a
|
||||
colour-scale legend is shown on the right. `plotters` writes the final
|
||||
PNG.
|
||||
4. **Render** — The base layer comes from Carto Voyager OSM tiles (fetched and
|
||||
cached separately for the basemap and for labels). Cloud cover is blended on top,
|
||||
then map labels are composited above the clouds. City markers, a title bar, and a
|
||||
colour legend complete the frame.
|
||||
|
||||
## Cloud cover colour scale
|
||||
|
||||
| Cloud cover | Colour |
|
||||
|---|---|
|
||||
| < 1 % | Land background (clear) |
|
||||
| 1–2 % | Very light blue |
|
||||
| 2–5 % | Light blue |
|
||||
| 5–10 % | Medium blue |
|
||||
| 10–20 % | Blue |
|
||||
| 20–50 % | Dark blue |
|
||||
| 50–100 % | Very dark blue |
|
||||
|
||||
## Marked locations
|
||||
|
||||
Falkensee · Pausin · Nauen · Hennigsdorf · Ketzin/Brückenkopf · Potsdam ·
|
||||
Berlin (Mitte)
|
||||
5. **Animate** — All frames are colour-quantised to 256-colour palettes (in parallel)
|
||||
and assembled into a looping GIF.
|
||||
|
||||
## Dependencies
|
||||
|
||||
All dependencies are pure Rust or vendored C code compiled into the binary.
|
||||
No system shared libraries are required at runtime beyond libc.
|
||||
All dependencies are pure Rust or vendored C compiled into the binary — no system
|
||||
shared libraries are required at runtime beyond libc.
|
||||
|
||||
| Crate | Role |
|
||||
|---|---|
|
||||
| `reqwest` (rustls-tls) | HTTP client with pure-Rust TLS |
|
||||
| `bzip2` | bzip2 decompression (vendors libbzip2) |
|
||||
| `plotters` + `plotters-bitmap` | PNG rendering with built-in bitmap font |
|
||||
| `bzip2` | Bzip2 decompression (vendors libbzip2) |
|
||||
| `plotters` + `plotters-bitmap` | PNG rendering |
|
||||
| `font8x8` | Built-in 8×8 bitmap font for map labels |
|
||||
| `clap` | CLI argument parsing |
|
||||
| `chrono` | Date/time handling |
|
||||
| `chrono` + `chrono-tz` | Date/time handling and timezone conversion |
|
||||
| `anyhow` | Error propagation |
|
||||
| `spade` | Delaunay triangulation for grid interpolation |
|
||||
| `rayon` | Parallel downloads, rendering, and GIF quantisation |
|
||||
| `image` | PNG decoding (loading frames for GIF assembly) |
|
||||
| `gif` | GIF encoding |
|
||||
|
||||
## Building
|
||||
|
||||
```sh
|
||||
cd cloud_cover_rs
|
||||
cargo build --release
|
||||
```
|
||||
|
||||
The compiled binary ends up at `target/release/cloud_cover`.
|
||||
The binary is at `target/release/cloud_cover`.
|
||||
|
||||
## Future ideas
|
||||
|
||||
- Configurable display timezone (currently hardcoded to Europe/Berlin)
|
||||
- Configurable city markers / observation sites
|
||||
- Configurable data sources
|
||||
- Vector map tile support for increased rendering resolution
|
||||
- Separation of cache and output files & automatic cache management
|
||||
- Automatic newest prediction fetching
|
||||
|
||||
BIN
docs/example.png
Normal file
BIN
docs/example.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 65 KiB |
52
src/grib.rs
52
src/grib.rs
@@ -171,20 +171,23 @@ fn decode_simple_packing(sec5: &[u8], sec7: &[u8], n: u32) -> Result<Vec<f32>> {
|
||||
|
||||
let r_bytes: [u8; 4] = sec5[0..4].try_into().unwrap();
|
||||
let r = f32::from_bits(u32::from_be_bytes(r_bytes));
|
||||
let e = read_i16(&sec5[4..6]);
|
||||
let d = read_i16(&sec5[6..8]);
|
||||
// GRIB2 encodes E and D as sign-and-magnitude (MSB = sign, lower 15 bits = magnitude),
|
||||
// NOT two's complement. Using i16::from_be_bytes would misinterpret e.g. 0x8009 as
|
||||
// -32759 instead of the correct -9, causing 2^E to underflow to zero.
|
||||
let e = read_grib_scale(&sec5[4..6]);
|
||||
let d = read_grib_scale(&sec5[6..8]);
|
||||
let n_bits = sec5[8] as u32;
|
||||
|
||||
if n_bits == 0 {
|
||||
// All values equal the reference value
|
||||
return Ok(vec![r * 10f32.powi(-(d as i32)); n as usize]);
|
||||
return Ok(vec![r * 10f32.powi(-d); n as usize]);
|
||||
}
|
||||
if n_bits > 32 {
|
||||
bail!("Implausible n_bits = {}", n_bits);
|
||||
}
|
||||
|
||||
let scale_e = 2f32.powi(e as i32);
|
||||
let scale_d = 10f32.powi(-(d as i32));
|
||||
let scale_e = 2f32.powi(e);
|
||||
let scale_d = 10f32.powi(-d);
|
||||
|
||||
let mut values = Vec::with_capacity(n as usize);
|
||||
let mut bit_offset: u64 = 0;
|
||||
@@ -239,8 +242,15 @@ fn read_u16(b: &[u8]) -> u16 {
|
||||
u16::from_be_bytes(b[0..2].try_into().unwrap())
|
||||
}
|
||||
|
||||
fn read_i16(b: &[u8]) -> i16 {
|
||||
i16::from_be_bytes(b[0..2].try_into().unwrap())
|
||||
/// Decode a GRIB2 sign-and-magnitude 16-bit integer.
|
||||
///
|
||||
/// The WMO GRIB2 specification encodes binary and decimal scale factors using
|
||||
/// a non-standard representation: bit 15 is the sign (1 = negative) and bits
|
||||
/// 14–0 hold the unsigned magnitude. This is different from two's complement.
|
||||
fn read_grib_scale(b: &[u8]) -> i32 {
|
||||
let raw = u16::from_be_bytes(b[0..2].try_into().unwrap());
|
||||
let magnitude = (raw & 0x7FFF) as i32;
|
||||
if raw & 0x8000 != 0 { -magnitude } else { magnitude }
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
@@ -270,6 +280,34 @@ mod tests {
|
||||
assert_eq!(result, vec![5.0f32, 5.0, 5.0]);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_read_grib_scale_sign_magnitude() {
|
||||
// 0x8009 = sign bit set, magnitude 9 → should be -9, not -32759 (two's complement)
|
||||
assert_eq!(read_grib_scale(&[0x80, 0x09]), -9);
|
||||
// 0x0009 = positive 9
|
||||
assert_eq!(read_grib_scale(&[0x00, 0x09]), 9);
|
||||
// 0x0000 = zero
|
||||
assert_eq!(read_grib_scale(&[0x00, 0x00]), 0);
|
||||
// 0x8000 = negative zero (sign bit set, magnitude 0) → 0
|
||||
assert_eq!(read_grib_scale(&[0x80, 0x00]), 0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_simple_packing_with_neg_e() {
|
||||
// Simulate DWD CLCT packing: R=0, E=-9 (0x8009), D=0, n_bits=16
|
||||
// A packed value of 51200 should decode to 51200 * 2^-9 = 100.0
|
||||
let r = 0.0f32;
|
||||
let mut sec5 = vec![0u8; 10];
|
||||
sec5[0..4].copy_from_slice(&r.to_bits().to_be_bytes());
|
||||
sec5[4..6].copy_from_slice(&[0x80, 0x09]); // E = -9 in sign-magnitude
|
||||
sec5[6..8].copy_from_slice(&[0x00, 0x00]); // D = 0
|
||||
sec5[8] = 16; // 16 bits per value
|
||||
// Packed value 51200 as big-endian u16
|
||||
let sec7 = [0xC8u8, 0x00]; // 0xC800 = 51200
|
||||
let result = decode_simple_packing(&sec5, &sec7, 1).unwrap();
|
||||
assert!((result[0] - 100.0f32).abs() < 0.01, "Expected ~100.0, got {}", result[0]);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_expand_with_bitmap() {
|
||||
// bitmap: 0b10110000 → positions 0, 2, 3 are present; 1, 4-7 are missing
|
||||
|
||||
173
src/main.rs
173
src/main.rs
@@ -2,10 +2,12 @@ mod download;
|
||||
mod grib;
|
||||
mod projection;
|
||||
mod render;
|
||||
mod tiles;
|
||||
|
||||
use anyhow::{bail, Context, Result};
|
||||
use chrono::{DateTime, Datelike, Timelike, Utc};
|
||||
use clap::Parser;
|
||||
use rayon::prelude::*;
|
||||
use std::path::PathBuf;
|
||||
|
||||
const BASE_URL: &str = "https://opendata.dwd.de/weather/nwp/icon-d2/grib/";
|
||||
@@ -28,6 +30,22 @@ 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,
|
||||
|
||||
/// Map centre latitude in degrees (default: 52.56)
|
||||
#[arg(long, value_name = "LAT", default_value_t = render::DEFAULT_CENTER_LAT)]
|
||||
center_lat: f64,
|
||||
|
||||
/// Map centre longitude in degrees (default: 13.08)
|
||||
#[arg(long, value_name = "LON", default_value_t = render::DEFAULT_CENTER_LON)]
|
||||
center_lon: f64,
|
||||
|
||||
/// OSM tile zoom level (default: 10)
|
||||
#[arg(long, value_name = "ZOOM", default_value_t = render::DEFAULT_OSM_ZOOM)]
|
||||
zoom: u8,
|
||||
}
|
||||
|
||||
fn main() -> Result<()> {
|
||||
@@ -77,26 +95,22 @@ fn main() -> Result<()> {
|
||||
dt.hour()
|
||||
);
|
||||
|
||||
// Download & decode coordinate grids
|
||||
// Build a dedicated thread pool for I/O-bound downloads, separate from
|
||||
// the global rayon pool used for CPU-bound rendering.
|
||||
let num_download_threads = (args.prediction_hours as usize + 2).min(16);
|
||||
let download_pool = rayon::ThreadPoolBuilder::new()
|
||||
.num_threads(num_download_threads)
|
||||
.build()
|
||||
.context("Failed to build download thread pool")?;
|
||||
|
||||
// Download coordinate grids (clat / clon) concurrently.
|
||||
eprintln!("Downloading coordinate grids (clat / clon)...");
|
||||
let lats = download_and_decode(
|
||||
&run_url,
|
||||
&date_str,
|
||||
"clat",
|
||||
"time-invariant",
|
||||
None,
|
||||
&cache_dir,
|
||||
)
|
||||
.context("clat")?;
|
||||
let lons = download_and_decode(
|
||||
&run_url,
|
||||
&date_str,
|
||||
"clon",
|
||||
"time-invariant",
|
||||
None,
|
||||
&cache_dir,
|
||||
)
|
||||
.context("clon")?;
|
||||
let (lats, lons) = download_pool.join(
|
||||
|| download_and_decode(&run_url, &date_str, "clat", "time-invariant", None, &cache_dir).context("clat"),
|
||||
|| download_and_decode(&run_url, &date_str, "clon", "time-invariant", None, &cache_dir).context("clon"),
|
||||
);
|
||||
let lats = lats?;
|
||||
let lons = lons?;
|
||||
|
||||
if lats.len() != lons.len() {
|
||||
bail!(
|
||||
@@ -107,21 +121,70 @@ fn main() -> Result<()> {
|
||||
}
|
||||
eprintln!(" Grid has {} points.", lats.len());
|
||||
|
||||
// Download, decode, and render each forecast hour
|
||||
// Download all forecast steps in parallel, then render sequentially.
|
||||
eprintln!("Downloading {} forecast frame(s) in parallel...", args.prediction_hours);
|
||||
let mut cloud_data: Vec<(u32, Vec<f32>)> = download_pool.install(|| {
|
||||
(0..args.prediction_hours)
|
||||
.into_par_iter()
|
||||
.map(|step| {
|
||||
let cloud = download_and_decode(
|
||||
&run_url,
|
||||
&date_str,
|
||||
"clct",
|
||||
"single-level",
|
||||
Some(step),
|
||||
&cache_dir,
|
||||
)
|
||||
.with_context(|| format!("clct step {}", step))?;
|
||||
Ok::<_, anyhow::Error>((step, cloud))
|
||||
})
|
||||
.collect::<Result<Vec<_>>>()
|
||||
})?;
|
||||
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 = render::make_projection(args.center_lat, args.center_lon, args.zoom);
|
||||
let tile_cache = tiles::fetch_tiles(
|
||||
&proj,
|
||||
render::MAP_W,
|
||||
render::MAP_H,
|
||||
args.zoom,
|
||||
&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 = render::make_projection(args.center_lat, args.center_lon, args.zoom);
|
||||
let label_cache = tiles::fetch_label_tiles(
|
||||
&proj,
|
||||
render::MAP_W,
|
||||
render::MAP_H,
|
||||
args.zoom,
|
||||
&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();
|
||||
|
||||
for step in 0..args.prediction_hours {
|
||||
let cloud = download_and_decode(
|
||||
&run_url,
|
||||
&date_str,
|
||||
"clct",
|
||||
"single-level",
|
||||
Some(step),
|
||||
&cache_dir,
|
||||
)
|
||||
.with_context(|| format!("clct step {}", step))?;
|
||||
|
||||
for (step, cloud) in cloud_data {
|
||||
if cloud.len() != lats.len() {
|
||||
bail!(
|
||||
"clct step {} has {} values but coordinate grids have {}",
|
||||
@@ -132,14 +195,15 @@ fn main() -> Result<()> {
|
||||
}
|
||||
|
||||
let forecast_dt = dt + chrono::Duration::hours(step as i64);
|
||||
let forecast_dt = forecast_dt.with_timezone(&chrono_tz::Europe::Berlin);
|
||||
let title = if step == 0 {
|
||||
format!(
|
||||
"Conditions at {} UTC",
|
||||
"Conditions at {} CET",
|
||||
forecast_dt.format("%Y-%m-%d %H:%M")
|
||||
)
|
||||
} else {
|
||||
format!(
|
||||
"Prediction at {} UTC (+{:02}h)",
|
||||
"Prediction at {} CET (+{:02}h)",
|
||||
forecast_dt.format("%Y-%m-%d %H:%M"),
|
||||
step
|
||||
)
|
||||
@@ -148,7 +212,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(), args.center_lat, args.center_lon, args.zoom)
|
||||
.with_context(|| format!("Rendering frame {}", step))?;
|
||||
|
||||
eprintln!(" [{:02}/{:02}] {}", step + 1, args.prediction_hours, out_name);
|
||||
@@ -160,6 +224,47 @@ fn main() -> Result<()> {
|
||||
eprintln!(" {}", p.display());
|
||||
}
|
||||
|
||||
if output_paths.len() > 1 {
|
||||
eprintln!("Building animation...");
|
||||
let gif_path = cache_dir.join("clct_animation.gif");
|
||||
create_gif(&output_paths, &gif_path).context("Creating GIF")?;
|
||||
eprintln!(" GIF: {}", gif_path.display());
|
||||
}
|
||||
|
||||
Ok(())
|
||||
}
|
||||
|
||||
fn create_gif(png_paths: &[PathBuf], out_path: &PathBuf) -> Result<()> {
|
||||
use gif::{Encoder, Frame, Repeat};
|
||||
|
||||
// Load PNGs and quantize to 256-colour palettes in parallel (the bottleneck).
|
||||
let n = png_paths.len();
|
||||
let frames: Vec<Frame<'static>> = png_paths
|
||||
.par_iter()
|
||||
.enumerate()
|
||||
.map(|(i, path)| -> Result<Frame<'static>> {
|
||||
let img = image::open(path)
|
||||
.with_context(|| format!("Loading PNG {}", path.display()))?
|
||||
.into_rgba8();
|
||||
let (width, height) = img.dimensions();
|
||||
let mut pixels = img.into_raw();
|
||||
// speed 8: slightly better quality than the default (10)
|
||||
let mut frame =
|
||||
Frame::from_rgba_speed(width as u16, height as u16, &mut pixels, 8);
|
||||
// Pause longer on first and last frame; 1 unit = 10 ms
|
||||
frame.delay = if i == 0 || i == n - 1 { 200 } else { 40 };
|
||||
Ok(frame)
|
||||
})
|
||||
.collect::<Result<Vec<_>>>()?;
|
||||
|
||||
// Write frames sequentially (GIF encoder is not thread-safe).
|
||||
let file = std::fs::File::create(out_path)?;
|
||||
let (w, h) = (frames[0].width, frames[0].height);
|
||||
let mut encoder = Encoder::new(file, w, h, &[])?;
|
||||
encoder.set_repeat(Repeat::Infinite)?;
|
||||
for frame in &frames {
|
||||
encoder.write_frame(frame)?;
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
|
||||
@@ -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.
|
||||
|
||||
345
src/render.rs
345
src/render.rs
@@ -2,31 +2,44 @@ 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;
|
||||
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;
|
||||
|
||||
/// Cloud-cover thresholds and their associated RGB colours.
|
||||
/// Values below the first threshold are drawn as the background colour.
|
||||
const LEVELS: [(f32, (u8, u8, u8)); 7] = [
|
||||
(1.0, (208, 232, 248)), // very thin
|
||||
(2.0, (160, 200, 232)),
|
||||
(5.0, (112, 168, 216)),
|
||||
(10.0, ( 64, 128, 192)),
|
||||
(20.0, ( 32, 80, 160)),
|
||||
(50.0, ( 16, 32, 96)),
|
||||
(100.0, ( 5, 10, 48)),
|
||||
];
|
||||
/// Default map image centre and OSM zoom level.
|
||||
pub const DEFAULT_CENTER_LAT: f64 = 52.56;
|
||||
pub const DEFAULT_CENTER_LON: f64 = 13.08;
|
||||
/// Default 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 DEFAULT_OSM_ZOOM: u8 = 10;
|
||||
|
||||
/// Build the orthographic projection whose viewport is sized so that the
|
||||
/// output resolution matches OSM tiles at `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(center_lat: f64, center_lon: f64, zoom: u8) -> OrthoProjection {
|
||||
let scale = f64::from(256u32 << 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(232, 236, 220); // muted green for land
|
||||
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);
|
||||
|
||||
@@ -41,16 +54,60 @@ const LOCATIONS: &[(&str, f64, f64)] = &[
|
||||
("Berlin (Mitte)", 13.39, 52.52),
|
||||
];
|
||||
|
||||
fn cloud_color(cover: f32) -> Option<RGBColor> {
|
||||
#[derive(Clone, Copy)]
|
||||
struct CloudVertex {
|
||||
position: Point2<f64>,
|
||||
cloud: f32,
|
||||
}
|
||||
|
||||
impl HasPosition for CloudVertex {
|
||||
type Scalar = f64;
|
||||
fn position(&self) -> Point2<f64> {
|
||||
self.position
|
||||
}
|
||||
}
|
||||
|
||||
/// Target colour for 100 % cloud cover: a light blue-white, like an overcast sky.
|
||||
/// Using a tinted target (rather than pure white) gives better contrast against
|
||||
/// the pale yellow-beige tones of the basemap.
|
||||
const CLOUD_TARGET: [f32; 3] = [226.0, 230.0, 239.0];
|
||||
|
||||
/// Boost the perceptual saturation of an RGB pixel by `factor`.
|
||||
/// Uses a luma-weighted approach: moves each channel away from the grey value.
|
||||
fn boost_saturation(pixel: [u8; 3], factor: f32) -> [u8; 3] {
|
||||
let [r, g, b] = pixel;
|
||||
let lum = 0.299 * r as f32 + 0.587 * g as f32 + 0.114 * b as f32;
|
||||
let clamp_u8 = |v: f32| v.round().clamp(0.0, 255.0) as u8;
|
||||
[
|
||||
clamp_u8(lum + factor * (r as f32 - lum)),
|
||||
clamp_u8(lum + factor * (g as f32 - lum)),
|
||||
clamp_u8(lum + factor * (b as f32 - lum)),
|
||||
]
|
||||
}
|
||||
|
||||
fn scale_cloud_cover(cover: f32) -> f32 {
|
||||
if cover < 1.0 {
|
||||
return None; // clear sky — use background
|
||||
0.0
|
||||
} else {
|
||||
30.0 + 0.7 * cover
|
||||
}
|
||||
for &(threshold, (r, g, b)) in LEVELS.iter().rev() {
|
||||
if cover >= threshold {
|
||||
return Some(RGBColor(r, g, b));
|
||||
}
|
||||
}
|
||||
|
||||
/// Map cloud cover (0–100 %) to a colour by blending the base color toward
|
||||
/// `CLOUD_TARGET`. Returns `None` for near-zero cover or NaN so the land background shows through.
|
||||
fn cloud_color(base_color: RGBColor, cover: f32) -> Option<RGBColor> {
|
||||
if cover.is_nan() || cover < 1.0 {
|
||||
return None;
|
||||
}
|
||||
Some(RGBColor(LEVELS[0].1 .0, LEVELS[0].1 .1, LEVELS[0].1 .2))
|
||||
let t = (cover / 100.0).clamp(0.0, 1.0);
|
||||
let blend = |base: u8, target: f32| -> u8 {
|
||||
(base as f32 + t * (target - base as f32)).round() as u8
|
||||
};
|
||||
Some(RGBColor(
|
||||
blend(base_color.0, CLOUD_TARGET[0]),
|
||||
blend(base_color.1, CLOUD_TARGET[1]),
|
||||
blend(base_color.2, CLOUD_TARGET[2]),
|
||||
))
|
||||
}
|
||||
|
||||
/// Draw ASCII text using the built-in 8×8 pixel font.
|
||||
@@ -90,6 +147,57 @@ fn draw_pixel_text(
|
||||
}
|
||||
}
|
||||
|
||||
/// 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.
|
||||
@@ -101,65 +209,151 @@ pub fn render_frame(
|
||||
lons: &[f32],
|
||||
cloud_cover: &[f32],
|
||||
title: &str,
|
||||
basemap: Option<&[[u8; 3]]>,
|
||||
labels: Option<&[[u8; 4]]>,
|
||||
center_lat: f64,
|
||||
center_lon: f64,
|
||||
zoom: u8,
|
||||
) -> Result<()> {
|
||||
let proj = OrthoProjection::new(CENTER_LAT, CENTER_LON, HALF_W_DEG, HALF_H_DEG);
|
||||
let proj = make_projection(center_lat, center_lon, zoom);
|
||||
|
||||
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 * 2;
|
||||
let map_w = MAP_W;
|
||||
let map_h = MAP_H;
|
||||
let n_pixels = map_w as usize * map_h as usize;
|
||||
|
||||
map_area.fill(&LAND_BG).context("fill map area")?;
|
||||
|
||||
// Plot grid points
|
||||
// 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 lat = lats[i] as f64;
|
||||
let lon = lons[i] as f64;
|
||||
|
||||
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 };
|
||||
|
||||
let cover = cloud_cover[i];
|
||||
if cover.is_nan() {
|
||||
continue; // missing value (bitmap) — keep background
|
||||
continue;
|
||||
}
|
||||
let color = match cloud_color(cover) {
|
||||
Some(c) => c,
|
||||
None => continue, // clear sky — keep land background
|
||||
};
|
||||
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 + 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
|
||||
|| row_f > map_h as f64 + TRI_MARGIN
|
||||
{
|
||||
continue;
|
||||
}
|
||||
let _ = tri.insert(CloudVertex { position: Point2::new(col_f, row_f), cloud: scale_cloud_cover(cover) });
|
||||
}
|
||||
|
||||
// Paint a 2×2 block to avoid gaps between grid points
|
||||
for dy in 0..2i32 {
|
||||
for dx in 0..2i32 {
|
||||
let c = col + dx;
|
||||
let r = row + dy;
|
||||
if c >= 0 && c < map_w as i32 && r >= 0 && r < map_h as i32 {
|
||||
map_area.draw_pixel((c, r), &color).ok();
|
||||
// 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);
|
||||
|
||||
// Build pixel buffer: start from basemap (saturation-boosted) or flat LAND_BG.
|
||||
const SAT_BOOST: f32 = 1.8;
|
||||
let mut pixel_buf: Vec<RGBColor> = if let Some(bm) = basemap {
|
||||
bm.iter().map(|&p| boost_saturation(p, SAT_BOOST)).map(|p| RGBColor(p[0], p[1], p[2])).collect()
|
||||
} else {
|
||||
vec![LAND_BG; n_pixels]
|
||||
};
|
||||
|
||||
// Blend cloud cover toward CLOUD_TARGET.
|
||||
for (idx, &cover) in cover_grid.iter().enumerate() {
|
||||
if let Some(cloud_col) = cloud_color(pixel_buf[idx], cover) {
|
||||
pixel_buf[idx] = cloud_col;
|
||||
}
|
||||
}
|
||||
|
||||
// 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();
|
||||
// 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 RGBColor(r, g, b) = pixel_buf[idx];
|
||||
pixel_buf[idx] = RGBColor(
|
||||
(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,
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
// Label offset slightly above and to the right of the marker
|
||||
draw_pixel_text(&map_area, label, col + 5, row - 9, 1, &BLACK);
|
||||
// Flush pixel buffer to drawing area.
|
||||
for row in 0..map_h as i32 {
|
||||
for col in 0..map_w as i32 {
|
||||
let pixel = pixel_buf[row as usize * map_w as usize + col as usize];
|
||||
map_area.draw_pixel((col, row), &pixel).ok();
|
||||
}
|
||||
}
|
||||
|
||||
// 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)
|
||||
@@ -204,22 +398,17 @@ fn draw_legend(
|
||||
);
|
||||
};
|
||||
|
||||
// "< 1%" row (clear sky)
|
||||
draw_entry("< 1%", LAND_BG, y as i32 + 22);
|
||||
|
||||
// Remaining threshold rows — use ASCII hyphen instead of en-dash
|
||||
let entries = [
|
||||
(1.0f32, "1-2%"),
|
||||
(2.0, "2-5%"),
|
||||
(5.0, "5-10%"),
|
||||
(10.0, "10-20%"),
|
||||
(20.0, "20-50%"),
|
||||
(50.0, "50-100%"),
|
||||
(0.0f32, "Clear"),
|
||||
(1.0, "1%"),
|
||||
(25.0, "25%"),
|
||||
(50.0, "50%"),
|
||||
(100.0, "100%"),
|
||||
];
|
||||
|
||||
for (i, &(lo, label)) in entries.iter().enumerate() {
|
||||
let row_y = y as i32 + 22 + (i as i32 + 1) * (box_h + 2);
|
||||
let color = cloud_color(lo).unwrap_or(BACKGROUND);
|
||||
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(LAND_BG, scale_cloud_cover(cover)).unwrap_or(LAND_BG);
|
||||
draw_entry(label, color, row_y);
|
||||
}
|
||||
|
||||
|
||||
354
src/tiles.rs
Normal file
354
src/tiles.rs
Normal file
@@ -0,0 +1,354 @@
|
||||
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)
|
||||
}
|
||||
|
||||
/// 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,
|
||||
zoom: u8,
|
||||
cache_dir: &Path,
|
||||
) -> Result<TileCache> {
|
||||
let (min_lon, max_lon, min_lat, max_lat) = map_bbox(proj, map_w, map_h);
|
||||
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,
|
||||
zoom: u8,
|
||||
cache_dir: &Path,
|
||||
) -> Result<LabelTileCache> {
|
||||
let (min_lon, max_lon, min_lat, max_lat) = map_bbox(proj, map_w, map_h);
|
||||
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