-6

configs/fetchWallpaper.nix

Reviewed

··· 1 1 - { 2 2 - lass.fetchWallpaper = { 3 3 - enable = true; 4 4 - url = "http://wallpaper.r/realwallpaper-krebs-stars-berlin.png"; 5 5 - }; 6 6 - }

-44

configs/realwallpaper.nix

Reviewed

··· 1 1 - { 2 2 - self, 3 3 - config, 4 4 - pkgs, 5 5 - ... 6 6 - }: 7 7 - { 8 8 - krebs.realwallpaper.enable = true; 9 9 - 10 10 - system.activationScripts.wallpaper-chmod = '' 11 11 - ${pkgs.coreutils}/bin/chmod +rx /var/realwallpaper 12 12 - ''; 13 13 - services.nginx.virtualHosts.wallpaper = { 14 14 - extraConfig = '' 15 15 - if ( $server_addr = "${ 16 16 - self.inputs.kartei.hosts.${config.networking.hostName}.nets.internet.ip4.addr 17 17 - }" ) { 18 18 - return 403; 19 19 - } 20 20 - ''; 21 21 - serverAliases = [ 22 22 - "wallpaper.r" 23 23 - ]; 24 24 - locations."/".extraConfig = '' 25 25 - autoindex on; 26 26 - root /var/realwallpaper/; 27 27 - ''; 28 28 - locations."/realwallpaper.png".extraConfig = '' 29 29 - root /var/realwallpaper/; 30 30 - ''; 31 31 - locations."/realwallpaper-krebs.png".extraConfig = '' 32 32 - root /var/realwallpaper/; 33 33 - ''; 34 34 - locations."/realwallpaper-krebs-stars.png".extraConfig = '' 35 35 - root /var/realwallpaper/; 36 36 - ''; 37 37 - locations."/realwallpaper-krebs-stars-berlin.png".extraConfig = '' 38 38 - root /var/realwallpaper/; 39 39 - ''; 40 40 - locations."/realwallpaper-video.mp4".extraConfig = '' 41 41 - root /var/realwallpaper/archive; 42 42 - ''; 43 43 - }; 44 44 - }

+1 -2

modules/default.nix

View file Reviewed

··· 4 4 ./fetchWallpaper.nix 5 5 ./sync-containers3.nix 6 6 ./parley.nix 7 7 - # copied from stockholm, see the header of each 7 7 + # copied from stockholm, see the header 8 8 ./reaktor2.nix 9 9 - ./realwallpaper.nix 10 9 ]; 11 10 }

-81

modules/realwallpaper.nix

Reviewed

··· 1 1 - # Copied from stockholm krebs/3modules/realwallpaper.nix (cc283503); 2 2 - # realwallpaper itself is pkgs/realwallpaper. 3 3 - { 4 4 - config, 5 5 - lib, 6 6 - pkgs, 7 7 - self, 8 8 - ... 9 9 - }: 10 10 - 11 11 - with lib; 12 12 - let 13 13 - cfg = config.krebs.realwallpaper; 14 14 - 15 15 - out = { 16 16 - options.krebs.realwallpaper = api; 17 17 - config = lib.mkIf cfg.enable imp; 18 18 - }; 19 19 - 20 20 - api = { 21 21 - enable = mkEnableOption "realwallpaper"; 22 22 - 23 23 - workingDir = mkOption { 24 24 - type = types.str; 25 25 - default = "/var/realwallpaper/"; 26 26 - }; 27 27 - 28 28 - marker = mkOption { 29 29 - type = types.str; 30 30 - default = "http://graph.r/network.json"; 31 31 - }; 32 32 - 33 33 - timerConfig = mkOption { 34 34 - type = types.unspecified; 35 35 - default = { 36 36 - OnCalendar = "*:0/15"; 37 37 - }; 38 38 - }; 39 39 - }; 40 40 - 41 41 - imp = { 42 42 - systemd.timers.realwallpaper = { 43 43 - description = "real wallpaper generator timer"; 44 44 - wantedBy = [ "timers.target" ]; 45 45 - 46 46 - timerConfig = cfg.timerConfig; 47 47 - }; 48 48 - 49 49 - systemd.services.realwallpaper = { 50 50 - description = "real wallpaper generator"; 51 51 - after = [ "network.target" ]; 52 52 - 53 53 - environment = { 54 54 - working_dir = cfg.workingDir; 55 55 - marker_url = cfg.marker; 56 56 - }; 57 57 - 58 58 - restartIfChanged = true; 59 59 - 60 60 - serviceConfig = { 61 61 - Type = "simple"; 62 62 - Restart = "on-failure"; 63 63 - ExecStart = "${ 64 64 - self.packages.${pkgs.stdenv.hostPlatform.system}.realwallpaper 65 65 - }/bin/generate-wallpaper"; 66 66 - User = "realwallpaper"; 67 67 - }; 68 68 - }; 69 69 - 70 70 - users.extraUsers.realwallpaper = { 71 71 - group = "realwallpaper"; 72 72 - home = cfg.workingDir; 73 73 - createHome = true; 74 74 - isSystemUser = true; 75 75 - }; 76 76 - 77 77 - users.groups.realwallpaper = { }; 78 78 - }; 79 79 - 80 80 - in 81 81 - out

-1

machines/mors/config.nix

View file Reviewed

··· 15 15 ../../configs/browsers.nix 16 16 ../../configs/pass.nix 17 17 ../../configs/steam.nix 18 18 - ../../configs/fetchWallpaper.nix 19 18 ../../configs/mail.nix 20 19 ../../configs/syncthing.nix 21 20 ../../configs/ableton.nix

-14

machines/neoprism/config.nix

View file Reviewed

··· 73 73 # where.is.lassul.us (location from home-assistant on a map) 74 74 ../../configs/websites/where.is.lassul.us 75 75 76 76 - # realwallpaper generator 77 77 - ../../configs/realwallpaper.nix 78 78 - 79 76 # paste + cyberlocker 80 77 ../../configs/paste 81 78 ··· 128 125 # lassul.us shouldn't be the default vhost here (nginx.nix already sets one) 129 126 services.nginx.virtualHosts."lassul.us" = { 130 127 default = lib.mkForce false; 131 131 - locations = { 132 132 - "= /wallpaper-marker.png".extraConfig = '' 133 133 - alias /var/realwallpaper/realwallpaper-marker.png; 134 134 - ''; 135 135 - "= /wallpaper.png".extraConfig = '' 136 136 - alias /var/realwallpaper/realwallpaper.png; 137 137 - ''; 138 138 - "= /wallpaper-stars-berlin.png".extraConfig = '' 139 139 - alias /var/realwallpaper/realwallpaper-krebs-stars-berlin.png; 140 140 - ''; 141 141 - }; 142 128 }; 143 129 144 130 # riot VM (HFOS): fresh libvirt guest on a dedicated bridged public IP.

-41

pkgs/realwallpaper/get_constellations.py

Reviewed

··· 1 1 - from astropy.coordinates import SkyCoord, ITRS, representation 2 2 - from astropy.time import Time 3 3 - import json 4 4 - import sys 5 5 - 6 6 - 7 7 - def convert_to_itrs(coord): 8 8 - c = SkyCoord(coord[0], coord[1], unit='degree', frame='icrs') 9 9 - c_itrs = c.transform_to(ITRS(obstime=Time.now())) 10 10 - rep = c_itrs.represent_as(representation.UnitSphericalRepresentation) 11 11 - return [rep.lat.deg, rep.lon.deg] 12 12 - 13 13 - 14 14 - def points_to_lines(points): 15 15 - lines = [] 16 16 - for x in range(len(points) - 1): 17 17 - lines.append([points[x], points[x+1]]) 18 18 - return lines 19 19 - 20 20 - 21 21 - def main(): 22 22 - with open(sys.argv[1]) as f: 23 23 - constellations = json.load(f)['features'] 24 24 - 25 25 - output = [] 26 26 - 27 27 - for const in constellations: 28 28 - for line in const['geometry']['coordinates']: 29 29 - transformed_line = [] 30 30 - for point in line: 31 31 - transformed_line.append(convert_to_itrs(point)) 32 32 - 33 33 - line_combined = points_to_lines(transformed_line) 34 34 - for l in line_combined: # noqa 35 35 - output.append(f'{l[0][0]} {l[0][1]} {l[1][0]} {l[1][1]} # {const["id"]}') # noqa 36 36 - 37 37 - print('\n'.join(output)) 38 38 - 39 39 - 40 40 - if __name__ == "__main__": 41 41 - main()

-24

pkgs/realwallpaper/grib2json.nix

Reviewed

··· 1 1 - # Copied from stockholm krebs/5pkgs/simple/grib2json (cc283503). 2 2 - { 3 3 - stdenv, 4 4 - fetchzip, 5 5 - jre, 6 6 - }: 7 7 - 8 8 - stdenv.mkDerivation { 9 9 - name = "grib2json"; 10 10 - src = fetchzip { 11 11 - url = "https://github.com/krebs/grib2json-bin/archive/31efe677b40fe491c988d50f96b59b1b7e2d46f7.zip"; 12 12 - hash = "sha256-yyijh5SQ0nMZqox8lweIwVKLHNl6aazmYWovcMTrfsA="; 13 13 - }; 14 14 - installPhase = '' 15 15 - mkdir -p $out/bin 16 16 - cp -r lib $out/ 17 17 - cat > $out/bin/grib2json << EOF 18 18 - #!/bin/sh 19 19 - set -x 20 20 - ${jre}/bin/java -jar $out/lib/grib2json-0.8.0-SNAPSHOT.jar "\$@" 21 21 - EOF 22 22 - chmod +x $out/bin/grib2json 23 23 - ''; 24 24 - }

-56

pkgs/realwallpaper/nomads-cloud.nix

Reviewed

··· 1 1 - # Copied from stockholm krebs/5pkgs/simple/nomads-cloud (cc283503). 2 2 - { 3 3 - callPackage, 4 4 - writers, 5 5 - coreutils, 6 6 - curl, 7 7 - jq, 8 8 - findutils, 9 9 - imagemagick, 10 10 - }: 11 11 - let 12 12 - grib2json = callPackage ./grib2json.nix { }; 13 13 - in 14 14 - writers.writeDashBin "nomads-cloud" '' 15 15 - prefix=$(mktemp -d) 16 16 - grib_path=$prefix/clouds.grib 17 17 - json_path=$prefix/clouds.json 18 18 - pgm_path=$prefix/clouds.pgm 19 19 - png_path=$1 20 20 - 21 21 - mkdir -p "$prefix" 22 22 - 23 23 - date=$(${coreutils}/bin/date +%Y%m%d) 24 24 - for hour in 18 12 06 00; do 25 25 - url="https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25_1hr.pl?file=gfs.t''${hour}z.pgrb2.0p25.anl&lev_entire_atmosphere_%5C%28considered_as_a_single_layer%5C%29=on&var_CWAT=on&leftlon=-180&rightlon=180&toplat=90&bottomlat=-90&dir=%2Fgfs.''${date}%2F''${hour}%2Fatmos" 26 26 - echo "$url" 27 27 - ${curl}/bin/curl -fsS "$url" > "$grib_path" 28 28 - if [ "$?" -eq 0 ]; then 29 29 - break 30 30 - fi 31 31 - done 32 32 - ${grib2json}/bin/grib2json --data "$grib_path" > "$json_path" 33 33 - 34 34 - width=$(${jq}/bin/jq '.[0].header.nx' < "$json_path") 35 35 - height=$(${jq}/bin/jq '.[0].header.ny' < "$json_path") 36 36 - 37 37 - # The maximum gray value. Must be bigger than 0 and less than 65536. 38 38 - maxval=1000 39 39 - 40 40 - # pgm - Netpbm grayscale image format 41 41 - # http://netpbm.sourceforge.net/doc/pgm.html 42 42 - { 43 43 - echo P2 44 44 - echo "$width $height" 45 45 - echo "$maxval" 46 46 - cat "$json_path" | 47 47 - ${jq}/bin/jq --argjson maxval "$maxval" -c ' 48 48 - ((.[0].data[]) * $maxval | round) 49 49 - ' | 50 50 - ${findutils}/bin/xargs -n "$width" 51 51 - } > "$pgm_path" 52 52 - 53 53 - ${imagemagick}/bin/convert -roll +50% "$pgm_path" "$png_path" 54 54 - 55 55 - rm -r "$prefix" 56 56 - ''

-374

pkgs/realwallpaper/package.nix

Reviewed

··· 1 1 - # Copied from stockholm krebs/5pkgs/simple/realwallpaper (cc283503). 2 2 - # nomads-cloud and grib2json are its stockholm dependencies, used only here. 3 3 - { pkgs, callPackage }: 4 4 - let 5 5 - nomads-cloud = callPackage ./nomads-cloud.nix { }; 6 6 - in 7 7 - (pkgs.writers.writeDashBin "generate-wallpaper" '' 8 8 - set -euf 9 9 - 10 10 - export PATH=${ 11 11 - pkgs.lib.makeBinPath ( 12 12 - [ 13 13 - nomads-cloud 14 14 - ] 15 15 - ++ (with pkgs; [ 16 16 - coreutils 17 17 - curl 18 18 - gnugrep 19 19 - gnused 20 20 - file 21 21 - findutils 22 22 - imagemagick 23 23 - inkscape 24 24 - jq 25 25 - xplanet 26 26 - ]) 27 27 - ) 28 28 - } 29 29 - 30 30 - # usage: getimg FILENAME URL 31 31 - fetch() { 32 32 - echo "fetch $1" 33 33 - curl -LsS -z "$1" -o "$1" "$2" 34 34 - } 35 35 - 36 36 - # usage: check_type FILENAME TYPE 37 37 - check_type() { 38 38 - if ! file -ib "$1" | grep -q "^$2/"; then 39 39 - echo "$1 is not of type $2" >&2 40 40 - rm "$1" 41 41 - return 1 42 42 - fi 43 43 - } 44 44 - 45 45 - # check if file exists and fetch only if missing 46 46 - fetch_once() { 47 47 - name=$1 48 48 - url=$2 49 49 - test -e "$name" || fetch "$name" "$url" 50 50 - } 51 51 - 52 52 - fetch_older_min() { 53 53 - min=$1 54 54 - name=$2 55 55 - url=$3 56 56 - if ! test "$(find $name -mmin -$min)"; then 57 57 - fetch "$name" "$url" 58 58 - fi 59 59 - } 60 60 - 61 61 - fetch_older_days() { 62 62 - days=$1 63 63 - name=$2 64 64 - url=$3 65 65 - if ! test "$(find $name -mtime -$days)"; then 66 66 - fetch "$name" "$url" 67 67 - fi 68 68 - } 69 69 - 70 70 - # usage: needs_rebuild DST SRC... 71 71 - needs_rebuild() { 72 72 - a="$1" 73 73 - shift 74 74 - if ! test -e "$a"; then 75 75 - #echo " $a does not exist" >&2 76 76 - result=0 77 77 - else 78 78 - result=1 79 79 - for b; do 80 80 - if check_type "$b" image; then 81 81 - if test "$b" -nt "$a"; then 82 82 - #echo " $b is newer than $a" >&2 83 83 - result=0 84 84 - fi 85 85 - fi 86 86 - done 87 87 - fi 88 88 - #case $result in 89 89 - # 0) echo "$a needs rebuild" >&2;; 90 90 - #esac 91 91 - return $result 92 92 - } 93 93 - 94 94 - get_neo_url() { 95 95 - url=$1 96 96 - curl -Ss "$url" | grep '3600 x 1800' | sed 's/.*href="\([^"]*\)".*/\1/' 97 97 - } 98 98 - 99 99 - main() { 100 100 - cd "''${working_dir:-$PWD}" 101 101 - 102 102 - # fetch source images in parallel 103 103 - fetch_once nightmap-raw.jpg \ 104 104 - 'https://eoimages.gsfc.nasa.gov/images/imagerecords/144000/144898/BlackMarble_2016_3km.jpg' & 105 105 - fetch_once daymap-raw.tif \ 106 106 - 'https://eoimages.gsfc.nasa.gov/images/imagerecords/57000/57752/land_shallow_topo_8192.tif' & 107 107 - 108 108 - fetch_once mercury-raw.svg \ 109 109 - 'https://upload.wikimedia.org/wikipedia/commons/2/2e/Mercury_symbol.svg' & 110 110 - fetch_once venus-raw.svg \ 111 111 - 'https://upload.wikimedia.org/wikipedia/commons/6/66/Venus_symbol.svg' & 112 112 - fetch_once mars-raw.svg \ 113 113 - 'https://upload.wikimedia.org/wikipedia/commons/b/b7/Mars_symbol.svg' & 114 114 - fetch_once jupiter-raw.svg \ 115 115 - 'https://upload.wikimedia.org/wikipedia/commons/2/26/Jupiter_symbol.svg' & 116 116 - fetch_once saturn-raw.svg \ 117 117 - 'https://upload.wikimedia.org/wikipedia/commons/7/74/Saturn_symbol.svg' & 118 118 - fetch_once uranus-raw.svg \ 119 119 - 'https://upload.wikimedia.org/wikipedia/commons/f/f1/Uranus_symbol.svg' & 120 120 - fetch_once neptune-raw.svg \ 121 121 - 'https://upload.wikimedia.org/wikipedia/commons/4/47/Neptune_symbol.svg' & 122 122 - 123 123 - fetch_once krebs-raw.svg \ 124 124 - 'https://raw.githubusercontent.com/krebs/painload/master/cholerab/bling/krebs_aquarium.svg' & 125 125 - 126 126 - fetch_older_min 720 ice-raw.jpg $(get_neo_url \ 127 127 - 'https://neo.gsfc.nasa.gov/view.php?datasetId=NISE_D') & 128 128 - fetch_older_days 1 snow-raw.jpg $(get_neo_url \ 129 129 - 'https://neo.gsfc.nasa.gov/view.php?datasetId=MOD10C1_E_SNOW') & 130 130 - fetch_older_days 1 chlora-raw.jpg $(get_neo_url \ 131 131 - 'https://neo.gsfc.nasa.gov/view.php?datasetId=MY1DMM_CHLORA') & 132 132 - fetch_older_days 1 fire-raw.jpg $(get_neo_url \ 133 133 - 'https://neo.gsfc.nasa.gov/view.php?datasetId=MOD14A1_E_FIRE') & 134 134 - 135 135 - # regular fetches 136 136 - fetch marker.json.tmp "''${marker_url:-}" || : 137 137 - if [ -s marker.json.tmp ]; then 138 138 - mv marker.json.tmp marker.json 139 139 - fi 140 140 - fetch sun-raw.jpg 'https://sdo.gsfc.nasa.gov/assets/img/latest/latest_512_0171.jpg' & 141 141 - 142 142 - wait 143 143 - 144 144 - # fetch clouds if they are older than 3h 145 145 - if ! test "$(find clouds-raw.png -mmin -180)"; then 146 146 - nomads-cloud clouds-raw.png 147 147 - fi 148 148 - 149 149 - in_size=3600x1800 150 150 - xplanet_out_size=3200x2500 151 151 - out_geometry=3200x1800+0+350 152 152 - 153 153 - for raw in \ 154 154 - nightmap-raw.jpg \ 155 155 - daymap-raw.tif \ 156 156 - chlora-raw.jpg \ 157 157 - clouds-raw.png \ 158 158 - ; 159 159 - do 160 160 - normal=''${raw%-raw.*}.png 161 161 - if needs_rebuild $normal $raw; then 162 162 - echo "make $normal; normalize $raw" >&2 163 163 - convert $raw -scale $in_size $normal 164 164 - fi 165 165 - done 166 166 - 167 167 - # remove snow from ice map 168 168 - if needs_rebuild ice.png ice-raw.jpg; then 169 169 - convert ice-raw.jpg -fuzz 20% -fill black -opaque white -colorspace gray -blur 0x6 -scale "$in_size" ice.png 170 170 - fi 171 171 - 172 172 - if needs_rebuild snow.png snow-raw.jpg; then 173 173 - convert snow-raw.jpg -fuzz 20% -fill '#DEDEDE' -opaque white -scale "$in_size" snow.png 174 174 - fi 175 175 - 176 176 - # make fire more red 177 177 - if needs_rebuild fire.png fire-raw.jpg; then 178 178 - convert fire-raw.jpg -fuzz 20% -fill '#ef840c' -opaque white -scale "$in_size" fire.png 179 179 - fi 180 180 - 181 181 - # cut out sun with alpha transparency 182 182 - if needs_rebuild sun.png sun-raw.jpg; then 183 183 - convert sun-raw.jpg \ 184 184 - \( +clone -colorspace HSB -fill white -draw "circle 256,256 256,54" -separate -delete 0,1 \) \ 185 185 - -compose copyopacity -composite -crop 512x472+0+20 -scale "100x100" sun.png 186 186 - fi 187 187 - 188 188 - if needs_rebuild krebs.png krebs-raw.svg; then 189 189 - inkscape --export-type="png" --export-width=16 --export-height=16 --export-filename=krebs.png krebs-raw.svg 190 190 - fi 191 191 - 192 192 - # -- Planets -- 193 193 - for planet in mercury venus mars jupiter saturn uranus neptune; do 194 194 - if needs_rebuild "$planet".png "$planet"-raw.svg; then 195 195 - sed -i 's/#000/#FE8019/g' "$planet"-raw.svg 196 196 - inkscape --export-type="png" --export-width=40 --export-height=40 --export-filename="$planet.png" "$planet-raw.svg" 197 197 - fi 198 198 - done 199 199 - 200 200 - # -- Daymap -- 201 201 - 202 202 - # merge with water chlora layer 203 203 - if needs_rebuild daymap-final.png daymap.png fire.png snow.png ice.png chlora.png; then 204 204 - convert daymap.png fire.png -compose lighten -composite daymap-1.png 205 205 - convert daymap-1.png ice.png -compose lighten -composite daymap-2.png 206 206 - convert daymap-2.png snow.png -compose lighten -composite daymap-3.png 207 207 - convert daymap-3.png chlora.png -compose lighten -composite daymap-final.png 208 208 - fi 209 209 - 210 210 - # -- Nightmap -- 211 211 - 212 212 - if needs_rebuild nightmap-final.png nightmap.png fire.png snow.png ice.png chlora.png; then 213 213 - convert nightmap.png fire.png -compose lighten -composite nightmap-1.png 214 214 - convert nightmap-1.png \( -fill black -colorize 70% ice.png \) -compose lighten -composite nightmap-2.png 215 215 - convert nightmap-2.png \( -fill black -colorize 70% snow.png \) -compose lighten -composite nightmap-3.png 216 216 - convert nightmap-3.png \( -fill black -colorize 70% chlora.png \) -compose lighten -composite nightmap-final.png 217 217 - fi 218 218 - 219 219 - # create marker file from json 220 220 - echo 'position=sun image=sun.png' > marker_file 221 221 - echo 'position=moon image=moon.png' >> marker_file 222 222 - echo 'position=mercury image=mercury.png' >> marker_file 223 223 - echo 'position=venus image=venus.png' >> marker_file 224 224 - echo 'position=mars image=mars.png' >> marker_file 225 225 - echo 'position=jupiter image=jupiter.png' >> marker_file 226 226 - echo 'position=saturn image=saturn.png' >> marker_file 227 227 - echo 'position=uranus image=uranus.png' >> marker_file 228 228 - echo 'position=neptune image=neptune.png' >> marker_file 229 229 - 230 230 - # generate moon 231 231 - xplanet -body moon --num_times 1 -origin earth \ 232 232 - --transpng moon.png --geometry 50x50 \ 233 233 - -config ${pkgs.writeText "moon.config" '' 234 234 - [moon] 235 235 - shade=10 236 236 - ''} 237 237 - 238 238 - # rebuild every time to update shadow 239 239 - xplanet --num_times 1 --geometry $xplanet_out_size \ 240 240 - --output xplanet-output.png --projection merc \ 241 241 - -config ${pkgs.writeText "xplanet.config" '' 242 242 - [earth] 243 243 - "Earth" 244 244 - map=daymap-final.png 245 245 - night_map=nightmap-final.png 246 246 - cloud_map=clouds.png 247 247 - cloud_threshold=1 248 248 - cloud_gamma=10 249 249 - shade=15 250 250 - ''} 251 251 - 252 252 - xplanet --num_times 1 --geometry $xplanet_out_size \ 253 253 - --output xplanet-marker-output.png --projection merc \ 254 254 - -config ${pkgs.writeText "xplanet-marker.config" '' 255 255 - [earth] 256 256 - "Earth" 257 257 - map=daymap-final.png 258 258 - night_map=nightmap-final.png 259 259 - cloud_map=clouds.png 260 260 - cloud_threshold=1 261 261 - cloud_gamma=10 262 262 - marker_file=marker_file 263 263 - shade=15 264 264 - ''} 265 265 - 266 266 - if [ -s marker.json ]; then 267 267 - jq -r 'to_entries[] | select(.value.latitude != null) | @json "\(.value.latitude) \(.value.longitude) image=krebs.png"' marker.json >> marker_file 268 268 - fi 269 269 - 270 270 - xplanet --num_times 1 --geometry $xplanet_out_size \ 271 271 - --output xplanet-krebs-output.png --projection merc \ 272 272 - -config ${pkgs.writeText "xplanet-krebs.config" '' 273 273 - [earth] 274 274 - "Earth" 275 275 - map=daymap-final.png 276 276 - night_map=nightmap-final.png 277 277 - cloud_map=clouds.png 278 278 - cloud_threshold=1 279 279 - cloud_gamma=10 280 280 - marker_file=marker_file 281 281 - shade=15 282 282 - ''} 283 283 - 284 284 - ${ 285 285 - pkgs.writers.writePython3 "get_constellations" { 286 286 - libraries = [ pkgs.python3Packages.astropy ]; 287 287 - } ./get_constellations.py 288 288 - } ${ 289 289 - pkgs.fetchurl { 290 290 - url = "https://raw.githubusercontent.com/ofrohn/d3-celestial/7e720a3de062059d4c5400a379146a601d9010e0/data/constellations.lines.json"; 291 291 - hash = "sha256-KU9mvvXVz1Cx4X8W0u+h2XoVExYSxo3ZNa3vbnNz4Tw="; 292 292 - } 293 293 - } > constellations.arcs || : # seems like astropy doesn't want to convert from icrs to itrs anymore 294 294 - 295 295 - xplanet --num_times 1 --geometry $xplanet_out_size \ 296 296 - --output xplanet-krebs-stars-output.png --projection merc \ 297 297 - -config ${pkgs.writeText "xplanet-krebs-stars.config" '' 298 298 - [default] 299 299 - 300 300 - arc_thickness=1 301 301 - arc_file=constellations.arcs 302 302 - 303 303 - [earth] 304 304 - "Earth" 305 305 - map=daymap-final.png 306 306 - night_map=nightmap-final.png 307 307 - cloud_map=clouds.png 308 308 - cloud_threshold=1 309 309 - cloud_gamma=10 310 310 - marker_file=marker_file 311 311 - shade=15 312 312 - ''} 313 313 - 314 314 - xplanet --num_times 1 --geometry $xplanet_out_size \ 315 315 - --latitude 52.520008 --longitude 13.404954 \ 316 316 - --output xplanet-krebs-stars-berlin-output.png --projection merc \ 317 317 - -config ${pkgs.writeText "xplanet-krebs-stars.config" '' 318 318 - [default] 319 319 - 320 320 - arc_thickness=1 321 321 - arc_file=constellations.arcs 322 322 - 323 323 - [earth] 324 324 - "Earth" 325 325 - map=daymap-final.png 326 326 - night_map=nightmap-final.png 327 327 - cloud_map=clouds.png 328 328 - cloud_threshold=1 329 329 - cloud_gamma=10 330 330 - marker_file=marker_file 331 331 - shade=15 332 332 - ''} 333 333 - 334 334 - # trim xplanet output 335 335 - if needs_rebuild realwallpaper.png xplanet-output.png; then 336 336 - convert xplanet-output.png -crop $out_geometry \ 337 337 - realwallpaper-tmp.png 338 338 - mv realwallpaper-tmp.png realwallpaper.png 339 339 - fi 340 340 - 341 341 - if needs_rebuild realwallpaper-marker.png xplanet-marker-output.png; then 342 342 - convert xplanet-marker-output.png -crop $out_geometry \ 343 343 - realwallpaper-marker-tmp.png 344 344 - mv realwallpaper-marker-tmp.png realwallpaper-marker.png 345 345 - fi 346 346 - 347 347 - if needs_rebuild realwallpaper-krebs.png xplanet-krebs-output.png; then 348 348 - convert xplanet-krebs-output.png -crop $out_geometry \ 349 349 - realwallpaper-krebs-tmp.png 350 350 - mv realwallpaper-krebs-tmp.png realwallpaper-krebs.png 351 351 - mkdir -p archive 352 352 - convert realwallpaper-krebs.png archive/"$(date -Is)".jpg 353 353 - fi 354 354 - 355 355 - if needs_rebuild realwallpaper-krebs-stars.png xplanet-krebs-stars-output.png; then 356 356 - convert xplanet-krebs-stars-output.png -crop $out_geometry \ 357 357 - realwallpaper-krebs-stars-tmp.png 358 358 - mv realwallpaper-krebs-stars-tmp.png realwallpaper-krebs-stars.png 359 359 - fi 360 360 - 361 361 - if needs_rebuild realwallpaper-krebs-stars-berlin.png xplanet-krebs-stars-berlin-output.png; then 362 362 - convert xplanet-krebs-stars-berlin-output.png -crop $out_geometry \ 363 363 - realwallpaper-krebs-stars-berlin-tmp.png 364 364 - mv realwallpaper-krebs-stars-berlin-tmp.png realwallpaper-krebs-stars-berlin.png 365 365 - fi 366 366 - } 367 367 - 368 368 - main "$@" 369 369 - '').overrideAttrs 370 370 - (old: { 371 371 - passthru = (old.passthru or { }) // { 372 372 - updateScript = ./update.sh; 373 373 - }; 374 374 - })

-65

pkgs/realwallpaper/update.sh

Reviewed

··· 1 1 - #!/usr/bin/env bash 2 2 - set -euo pipefail 3 3 - 4 4 - # realwallpaper has no version of its own; it pins two upstream GitHub assets 5 5 - # by commit. Bump each to the head of its default branch and refresh its hash. 6 6 - # Run by update-pkgs from the flake root ($0 is a read-only store copy). 7 7 - attr=realwallpaper 8 8 - dir=pkgs/realwallpaper 9 9 - system=$(nix eval --raw --impure --expr 'builtins.currentSystem') 10 10 - build_attr=".#packages.$system.$attr" 11 11 - 12 12 - # file | owner/repo | branch | URL prefix up to (and excluding) the 40-char rev 13 13 - assets=( 14 14 - "$dir/grib2json.nix|krebs/grib2json-bin|master|https://github.com/krebs/grib2json-bin/archive/" 15 15 - "$dir/package.nix|ofrohn/d3-celestial|master|https://raw.githubusercontent.com/ofrohn/d3-celestial/" 16 16 - ) 17 17 - 18 18 - # The probe build is *expected* to fail (hash mismatch); `|| true` keeps 19 19 - # pipefail/errexit from turning that into an empty `got`. awk must not 20 20 - # `exit` early or tee dies of SIGPIPE and takes the pipeline with it. 21 21 - probe_hash() { 22 22 - { nix build --no-link "$build_attr" 2>&1 || true; } | tee /dev/stderr \ 23 23 - | awk '/got:/ && !seen { print $2; seen = 1 }' || true 24 24 - } 25 25 - 26 26 - changed=0 27 27 - for asset in "${assets[@]}"; do 28 28 - IFS='|' read -r file repo branch prefix <<<"$asset" 29 29 - prefix_re=${prefix//./\\.} 30 30 - 31 31 - old_rev=$(grep -oE "${prefix_re}[0-9a-f]{40}" "$file" | sed "s|^${prefix_re}||") 32 32 - new_rev=$(curl -fsSL "https://api.github.com/repos/$repo/commits/$branch" | jq -r .sha) 33 33 - if [ -z "$old_rev" ] || [ -z "$new_rev" ] || [ "$new_rev" = null ]; then 34 34 - echo "$attr: could not determine revisions for $repo (old='$old_rev' new='$new_rev')" >&2 35 35 - exit 1 36 36 - fi 37 37 - if [ "$old_rev" = "$new_rev" ]; then 38 38 - echo "$attr: $repo already at $new_rev" 39 39 - continue 40 40 - fi 41 41 - 42 42 - # Bump the rev and blank the hash that follows the URL: with the old hash in 43 43 - # place nix would reuse the cached FOD and "succeed" on stale content. 44 44 - sed -i \ 45 45 - -e "s|${prefix_re}${old_rev}|${prefix}${new_rev}|" \ 46 46 - -e "\|${prefix_re}${new_rev}|,\|hash = |s|hash = \"[^\"]*\"|hash = \"\"|" \ 47 47 - "$file" 48 48 - 49 49 - # Only this asset's hash is blank, so the first mismatch belongs to it. 50 50 - got=$(probe_hash) 51 51 - if [ -z "$got" ]; then 52 52 - echo "$attr: no hash mismatch reported for $repo" >&2 53 53 - exit 1 54 54 - fi 55 55 - sed -i "\|${prefix_re}${new_rev}|,\|hash = |s|hash = \"\"|hash = \"$got\"|" "$file" 56 56 - echo "$attr: $repo $old_rev → $new_rev ($got)" 57 57 - changed=1 58 58 - done 59 59 - 60 60 - if [ "$changed" = 0 ]; then 61 61 - echo "$attr: already up to date" 62 62 - exit 0 63 63 - fi 64 64 - 65 65 - nix build --no-link "$build_attr"

+43

tools/realwallpaper/flake-module.nix

View file Reviewed

··· 1 1 + { 2 2 + perSystem = 3 3 + { pkgs, ... }: 4 4 + let 5 5 + python = pkgs.python3.withPackages ( 6 6 + ps: with ps; [ 7 7 + astropy 8 8 + numpy 9 9 + pillow 10 10 + ] 11 11 + ); 12 12 + constellations = pkgs.fetchurl { 13 13 + url = "https://raw.githubusercontent.com/ofrohn/d3-celestial/7e720a3de062059d4c5400a379146a601d9010e0/data/constellations.lines.json"; 14 14 + hash = "sha256-KU9mvvXVz1Cx4X8W0u+h2XoVExYSxo3ZNa3vbnNz4Tw="; 15 15 + }; 16 16 + krebsIcon = 17 17 + pkgs.runCommand "krebs.png" 18 18 + { 19 19 + nativeBuildInputs = [ pkgs.librsvg ]; 20 20 + src = pkgs.fetchurl { 21 21 + url = "https://raw.githubusercontent.com/krebs/painload/25bb340f4d8d05bb3ab6792277951b77923e7010/cholerab/bling/krebs_aquarium.svg"; 22 22 + hash = "sha256-KLoY58M3bPcI0WrqomPbDNewKz+akeR0oujccDeUqnY="; 23 23 + }; 24 24 + } 25 25 + '' 26 26 + rsvg-convert -w 64 -h 64 "$src" -o "$out" 27 27 + ''; 28 28 + in 29 29 + { 30 30 + packages.realwallpaper = 31 31 + (pkgs.writeShellApplication { 32 32 + name = "realwallpaper"; 33 33 + runtimeInputs = [ python ]; 34 34 + text = '' 35 35 + export REALWALLPAPER_CONSTELLATIONS=${constellations} 36 36 + export REALWALLPAPER_KREBS_ICON=${krebsIcon} 37 37 + export REALWALLPAPER_FONT=${pkgs.dejavu_fonts}/share/fonts/truetype/DejaVuSans.ttf 38 38 + exec python3 ${./realwallpaper.py} "$@" 39 39 + ''; 40 40 + }).overrideAttrs 41 41 + { passthru.usage = builtins.readFile ./usage.kdl; }; 42 42 + }; 43 43 + }

+1386

tools/realwallpaper/realwallpaper.py

View file Reviewed

··· 1 1 + """Render a live world map wallpaper: oblique Mercator with a lens on home. 2 2 + 3 3 + The globe is rotated so home (default Berlin) is in the middle of the screen 4 4 + on the projection's equator, then shown in Mercator filling the screen edge 5 5 + to edge. Around home a smooth radial lens magnifies the map and blends back 6 6 + into the plain map at its rim. 7 7 + 8 8 + Imagery is fetched over WMS at the resolution each pixel needs: a few local 9 9 + levels around home, one global level, and downsampled global mipmaps. The 10 10 + night side is the day map darkened, plus city lights as a glow. Every fetch 11 11 + is cached with a TTL; a failing source falls back to its last good copy or 12 12 + is skipped, so one dead upstream no longer kills the whole wallpaper. 13 13 + """ 14 14 + 15 15 + from __future__ import annotations 16 16 + 17 17 + import argparse 18 18 + import dataclasses 19 19 + import concurrent.futures 20 20 + import datetime as dt 21 21 + import hashlib 22 22 + import io 23 23 + import json 24 24 + import math 25 25 + import os 26 26 + import sys 27 27 + import time 28 28 + import urllib.parse 29 29 + import urllib.request 30 30 + import warnings 31 31 + from dataclasses import dataclass, field 32 32 + from pathlib import Path 33 33 + from typing import Callable 34 34 + 35 35 + import numpy as np 36 36 + from PIL import Image, ImageDraw, ImageFilter, ImageFont 37 37 + 38 38 + R_EARTH = 6371.0088 # km 39 39 + KM_PER_DEG = math.pi * R_EARTH / 180 40 40 + GLOBAL_WIDTH = 8192 # px of the global equirectangular level 41 41 + # px per WMS request side. The servers render on the fly at a roughly fixed 42 42 + # px/s per request (EOX: ~0.35 Mpx/s) but serve requests in parallel, so 43 43 + # many medium tiles fetched concurrently are much faster than few big ones 44 44 + MAX_TILE = 2048 45 45 + TILE_WORKERS = 4 # concurrent requests per level 46 46 + MAX_LEVEL_SIDE = 8192 # px per local level side 47 47 + USER_AGENT = "realwallpaper (+https://github.com/Lassulus/superconfig)" 48 48 + 49 49 + GIBS = "https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi" 50 50 + EUMETSAT = "https://view.eumetsat.int/geoserver/wms" 51 51 + EOX = "https://tiles.maps.eox.at/wms" 52 52 + 53 53 + 54 54 + def log(msg: str) -> None: 55 55 + print(msg, file=sys.stderr, flush=True) 56 56 + 57 57 + 58 58 + def smoothstep(e0: float, e1: float, x: np.ndarray) -> np.ndarray: 59 59 + t = np.clip((x - e0) / (e1 - e0), 0.0, 1.0) 60 60 + return t * t * (3.0 - 2.0 * t) 61 61 + 62 62 + 63 63 + # -------------------------------------------------------------------------- 64 64 + # projection 65 65 + 66 66 + 67 67 + def solve_cylinder(half_height: float) -> float: 68 68 + """Parameter a of the cylindrical y(lat) = ln tan(pi/4 + a lat/2) / a 69 69 + whose poles sit exactly half_height (radians of longitude) from the 70 70 + equator. a -> 0 is plate carrée (pi/2), a = 1 is Mercator (infinite); 71 71 + Miller is a = 0.8. Wider screens than 2:1 fall back to plate carrée.""" 72 72 + lo, hi = 1e-6, 1 - 1e-9 73 73 + if half_height <= math.pi / 2: 74 74 + return lo 75 75 + for _ in range(200): 76 76 + a = (lo + hi) / 2 77 77 + if math.log(math.tan(math.pi / 4 + a * math.pi / 4)) / a > half_height: 78 78 + hi = a 79 79 + else: 80 80 + lo = a 81 81 + return (lo + hi) / 2 82 82 + 83 83 + 84 84 + def rotate_to_home(lat0: float, lat, dlon): 85 85 + """Geographic (lat, lon - lon0) -> frame where home is lat' = lon' = 0. 86 86 + 87 87 + Rotating about the east-west axis keeps local north pointing up at home. 88 88 + All angles in radians. 89 89 + """ 90 90 + x = np.cos(lat) * np.cos(dlon) 91 91 + y = np.cos(lat) * np.sin(dlon) 92 92 + z = np.sin(lat) 93 93 + s0, c0 = math.sin(lat0), math.cos(lat0) 94 94 + xr = x * c0 + z * s0 95 95 + zr = -x * s0 + z * c0 96 96 + return np.arcsin(np.clip(zr, -1.0, 1.0)), np.arctan2(y, xr) 97 97 + 98 98 + 99 99 + def rotate_from_home(lat0: float, lat_r, lon_r): 100 100 + """Inverse of rotate_to_home: returns (lat, lon - lon0) in radians.""" 101 101 + xr = np.cos(lat_r) * np.cos(lon_r) 102 102 + y = np.cos(lat_r) * np.sin(lon_r) 103 103 + zr = np.sin(lat_r) 104 104 + s0, c0 = math.sin(lat0), math.cos(lat0) 105 105 + x = xr * c0 - zr * s0 106 106 + z = xr * s0 + zr * c0 107 107 + return np.arcsin(np.clip(z, -1.0, 1.0)), np.arctan2(y, x) 108 108 + 109 109 + 110 110 + class LensMap: 111 111 + """Screen-filling oblique cylindrical map around home, with a lens. 112 112 + 113 113 + The globe is rotated so home sits on the projection's equator in the 114 114 + middle of the screen (local north up); its poles land top and bottom 115 115 + only when home is on the real equator. The cylinder is Mercator-like 116 116 + (conformal along that equator) but with finite poles, y = ln tan(pi/4 + 117 117 + a lat/2) / a, with a picked so the whole globe fits the screen exactly: 118 118 + nothing is cut off, only squeezed towards the top and bottom edges. 119 119 + 120 120 + Projections (--projection): 121 121 + - lens: the above plus the magnifying lens around home (default) 122 122 + - plain: the above without the lens, no distortion beyond the cylinder 123 123 + - mercator: the classic map, north up, true Mercator centred on the 124 124 + equator and the home longitude; the poles are cut off (about +-71 125 125 + degrees fit 16:9) and there is no lens 126 126 + 127 127 + The lens is focus + glue + context (Yamamoto et al. 2009; the same idea 128 128 + as Snyder's 1987 "magnifying-glass" azimuthal projections): a lens 129 129 + "radius" n is measured as a circle, or following the screen's aspect 130 130 + as an ellipse (p = 2) or rounded rectangle (p = 4), the norm 131 131 + 132 132 + n = ((|dx| / ax)^p + (|dy| / ay)^p)^(1/p) 133 133 + 134 134 + (ax, ay = half the screen size times --lens). A screen point at n shows 135 135 + the map point at h(n) = n / z(n) along the same ray, where z is the 136 136 + local zoom: M in the focus (n < f, undistorted), 1 from n = 1 on (plain 137 137 + map), and in between, with t = (n - f) / (1 - f), 138 138 + 139 139 + ln z(t) = ln M * (1 - I(t) / I(1)), I(t) = int_0^t u^(a-1) (1 - u) du 140 140 + 141 141 + How squashed shapes are is the ratio of radial to tangential scale, 142 142 + 1 - t z'/z; its excess over 1 is proportional to t^a (1 - t): none at home, largest at 143 143 + t = a / (a + 1) =: --falloff, back to 1 at the lens edge. So --falloff 144 144 + says where the distortion sits: low keeps it near home, high pushes it 145 145 + out towards the border and leaves the middle closer to a plain zoom. 146 146 + z never increases outwards, so h' > 0 and the lens never folds. 147 147 + """ 148 148 + 149 149 + def __init__( 150 150 + self, 151 151 + lat: float, 152 152 + lon: float, 153 153 + width: int, 154 154 + height: int, 155 155 + magnification: float, 156 156 + lens: float, 157 157 + shape: str, 158 158 + focus: float, 159 159 + falloff: float, 160 160 + projection: str = "lens", 161 161 + ) -> None: 162 162 + if magnification < 1: 163 163 + raise SystemExit("--magnification must be >= 1") 164 164 + if not 0 <= focus < 1: 165 165 + raise SystemExit("--focus must be in [0, 1)") 166 166 + if not 0 < falloff < 1: 167 167 + raise SystemExit("--falloff must be in (0, 1)") 168 168 + self.lat0 = math.radians(lat) 169 169 + self.lon0 = math.radians(lon) 170 170 + self.width = width 171 171 + self.height = height 172 172 + self.m = magnification if projection == "lens" else 1.0 173 173 + self.k = 2 * math.pi / width # radians of longitude per px 174 174 + if projection == "mercator": 175 175 + self.rot_lat = 0.0 # no rotation: the real equator is the axis 176 176 + self.a = 1.0 # true Mercator 177 177 + else: 178 178 + self.rot_lat = self.lat0 # home on the projection's equator 179 179 + self.a = solve_cylinder(math.pi * height / width) 180 180 + self.px = width / 2 181 181 + self.py = height / 2 182 182 + self.shape = shape 183 183 + if shape == "circle": 184 184 + self.ax = self.ay = lens * height / 2 185 185 + else: 186 186 + self.ax = lens * width / 2 187 187 + self.ay = lens * height / 2 188 188 + # tabulate the profile; forward() inverts it by interpolation 189 189 + self._s = np.linspace(0.0, 1.0, 8001) 190 190 + t = np.clip((self._s - focus) / (1 - focus), 0.0, 1.0) 191 191 + a = falloff / (1 - falloff) 192 192 + with np.errstate(divide="ignore", invalid="ignore"): 193 193 + f_int = np.where(t > 0, t ** (a - 1) * (1 - t), 0.0) 194 194 + integral = np.concatenate([[0.0], np.cumsum((f_int[1:] + f_int[:-1]) / 2 * np.diff(t))]) 195 195 + zoom = np.exp(math.log(self.m) * (1 - integral / max(integral[-1], 1e-12))) 196 196 + self._h = self._s / zoom 197 197 + self._dh = np.gradient(self._h, self._s) 198 198 + # sample the screen to plan which resolutions are needed where 199 199 + ys, xs = np.mgrid[0:height:6, 0:width:6] + 3.0 200 200 + lat_g, lon_g, q = self.inverse(xs.ravel(), ys.ravel()) 201 201 + s = self.norm(xs.ravel() - self.px, ys.ravel() - self.py) 202 202 + self._plan_q = q[s < 1] 203 203 + self._plan_d = great_circle_km(lat, lon, lat_g[s < 1], lon_g[s < 1]) 204 204 + self.q0 = float(q.min()) 205 205 + self.q_max = float(q.max()) 206 206 + 207 207 + def norm(self, dx, dy): 208 208 + if self.shape == "circle": 209 209 + return np.hypot(dx, dy) / self.ax 210 210 + if self.shape == "ellipse": 211 211 + return np.hypot(dx / self.ax, dy / self.ay) 212 212 + return ((np.abs(dx) / self.ax) ** 4 + (np.abs(dy) / self.ay) ** 4) ** 0.25 213 213 + 214 214 + def h(self, s): 215 215 + return np.interp(s, self._s, self._h) 216 216 + 217 217 + def dh(self, s): 218 218 + return np.interp(s, self._s, self._dh) 219 219 + 220 220 + def distance_for_demand(self, q: float) -> float: 221 221 + """Radius (km) around home inside which some pixel needs finer than q.""" 222 222 + need = self._plan_d[self._plan_q < q] 223 223 + return float(min(need.max(), 3500.0)) if len(need) else 0.0 224 224 + 225 225 + def inverse(self, xs: np.ndarray, ys: np.ndarray): 226 226 + """Screen px -> lat, lon (degrees) and ground size of the px (km).""" 227 227 + dx, dy = xs - self.px, ys - self.py 228 228 + s = self.norm(dx, dy) 229 229 + inside = s < 1 230 230 + si = np.where(inside, s, 1.0) 231 231 + with np.errstate(divide="ignore", invalid="ignore"): 232 232 + tangential = np.where(si > 1e-9, self.h(si) / si, 1 / self.m) 233 233 + radial = self.dh(si) 234 234 + scale = np.where(inside, tangential, 1.0) 235 235 + wx = self.px + dx * scale 236 236 + wy = self.py + dy * scale 237 237 + lon_r = (wx - self.px) * self.k 238 238 + u = (self.py - wy) * self.k 239 239 + lat_r = (2 / self.a) * (np.arctan(np.exp(self.a * u)) - math.pi / 4) 240 240 + lat_r = np.clip(lat_r, -math.pi / 2, math.pi / 2) 241 241 + lat, dlon = rotate_from_home(self.rot_lat, lat_r, lon_r) 242 242 + lon = (self.lon0 + dlon + math.pi) % (2 * math.pi) - math.pi 243 243 + # the lens squeezes radially more than tangentially in its ring; use 244 244 + # the finer scale for sharpness, but not below a quarter of the other 245 245 + fine = np.minimum(tangential, radial) 246 246 + coarse = np.maximum(tangential, radial) 247 247 + lens_f = np.where(inside, np.maximum(fine, coarse / 4), 1.0) 248 248 + # cylinder: east-west px size R k cos(lat), north-south R k cos(a lat) 249 249 + base_f = np.maximum(np.cos(lat_r), np.cos(self.a * lat_r) / 4) 250 250 + q = R_EARTH * self.k * base_f * lens_f 251 251 + return np.degrees(lat), np.degrees(lon), np.maximum(q, 1e-3) 252 252 + 253 253 + def forward(self, lat_deg, lon_deg): 254 254 + lat = np.radians(np.asarray(lat_deg, dtype=np.float64)) 255 255 + dlon = np.radians(np.asarray(lon_deg, dtype=np.float64)) - self.lon0 256 256 + lat_r, lon_r = rotate_to_home(self.rot_lat, lat, dlon) 257 257 + lat_r = np.clip(lat_r, -math.radians(89.5), math.radians(89.5)) 258 258 + wx = self.px + lon_r / self.k 259 259 + wy = self.py - np.log(np.tan(math.pi / 4 + self.a * lat_r / 2)) / (self.a * self.k) 260 260 + dx, dy = wx - self.px, wy - self.py 261 261 + r = self.norm(dx, dy) 262 262 + s = np.interp(np.minimum(r, 1.0), self._h, self._s) 263 263 + with np.errstate(divide="ignore", invalid="ignore"): 264 264 + scale = np.where(r < 1, np.where(r > 1e-12, s / r, self.m), 1.0) 265 265 + return self.px + dx * scale, self.py + dy * scale 266 266 + 267 267 + 268 268 + def great_circle_km(lat0: float, lon0: float, lat, lon) -> np.ndarray: 269 269 + la0, lo0 = math.radians(lat0), math.radians(lon0) 270 270 + la, lo = np.radians(lat), np.radians(lon) 271 271 + c = math.sin(la0) * np.sin(la) + math.cos(la0) * np.cos(la) * np.cos(lo - lo0) 272 272 + return np.arccos(np.clip(c, -1.0, 1.0)) * R_EARTH 273 273 + 274 274 + 275 275 + # -------------------------------------------------------------------------- 276 276 + # sources and levels 277 277 + 278 278 + 279 279 + @dataclass(frozen=True) 280 280 + class Wms: 281 281 + key: str 282 282 + base: str 283 283 + layer: str 284 284 + native_km: float 285 285 + ttl: float 286 286 + fmt: str = "image/png" 287 287 + time: str | None = None 288 288 + # TIME=<today - n days> for each n, newest first; gaps are filled from 289 289 + # older days (daily polar-orbiter products are full of swath gaps) 290 290 + days_back: tuple[int, ...] = () 291 291 + # fetch this many times finer and box-filter down: some servers (EOX) 292 292 + # pick single pixels when asked for a coarse scale, which aliases 293 293 + oversample: int = 1 294 294 + # per-channel gain to match the colours of the other day source 295 295 + gain: tuple[float, float, float] | None = None 296 296 + 297 297 + def times(self, when: dt.datetime) -> list[str | None]: 298 298 + if self.days_back: 299 299 + return [(when - dt.timedelta(days=n)).strftime("%Y-%m-%d") for n in self.days_back] 300 300 + return [self.time] 301 301 + 302 302 + def url(self, bbox: tuple[float, float, float, float], w: int, h: int, t: str | None) -> str: 303 303 + s, west, n, e = bbox 304 304 + params = { 305 305 + "SERVICE": "WMS", 306 306 + "VERSION": "1.3.0", 307 307 + "REQUEST": "GetMap", 308 308 + "LAYERS": self.layer, 309 309 + "STYLES": "", 310 310 + "CRS": "EPSG:4326", 311 311 + "BBOX": f"{s:.6f},{west:.6f},{n:.6f},{e:.6f}", 312 312 + "WIDTH": str(w), 313 313 + "HEIGHT": str(h), 314 314 + "FORMAT": self.fmt, 315 315 + "TRANSPARENT": "TRUE", 316 316 + } 317 317 + if t: 318 318 + params["TIME"] = t 319 319 + return self.base + "?" + urllib.parse.urlencode(params) 320 320 + 321 321 + 322 322 + @dataclass 323 323 + class Level: 324 324 + res: float # km per px 325 325 + south: float 326 326 + west: float 327 327 + north: float 328 328 + east: float 329 329 + width: int 330 330 + height: int 331 331 + wrap: bool 332 332 + data: np.ndarray | None = None 333 333 + 334 334 + @property 335 335 + def bbox(self) -> tuple[float, float, float, float]: 336 336 + return (self.south, self.west, self.north, self.east) 337 337 + 338 338 + 339 339 + def local_level(proj: LensMap, res: float, radius_km: float) -> Level: 340 340 + lat = math.degrees(proj.lat0) 341 341 + lon = math.degrees(proj.lon0) 342 342 + dlat = math.degrees(radius_km / R_EARTH) 343 343 + south, north = max(-90.0, lat - dlat), min(90.0, lat + dlat) 344 344 + coslat = math.cos(proj.lat0) 345 345 + s = math.sin(radius_km / R_EARTH) / coslat if coslat > 1e-6 else 2.0 346 346 + if north >= 90.0 or south <= -90.0 or s >= 1.0: 347 347 + west, east = -180.0, 180.0 348 348 + else: 349 349 + dlon = math.degrees(math.asin(s)) 350 350 + west, east = lon - dlon, lon + dlon 351 351 + mid = math.cos(math.radians((south + north) / 2)) 352 352 + width_km = (east - west) * KM_PER_DEG * max(mid, 0.05) 353 353 + height_km = (north - south) * KM_PER_DEG 354 354 + res = max(res, width_km / MAX_LEVEL_SIDE, height_km / MAX_LEVEL_SIDE) 355 355 + return Level( 356 356 + res=res, 357 357 + south=south, 358 358 + west=west, 359 359 + north=north, 360 360 + east=east, 361 361 + width=max(2, math.ceil(width_km / res)), 362 362 + height=max(2, math.ceil(height_km / res)), 363 363 + wrap=False, 364 364 + ) 365 365 + 366 366 + 367 367 + def global_level(width: int) -> Level: 368 368 + return Level( 369 369 + res=360 / width * KM_PER_DEG, 370 370 + south=-90.0, 371 371 + west=-180.0, 372 372 + north=90.0, 373 373 + east=180.0, 374 374 + width=width, 375 375 + height=width // 2, 376 376 + wrap=True, 377 377 + ) 378 378 + 379 379 + 380 380 + def plan_local(proj: LensMap, native_km: float) -> list[Level]: 381 381 + """Local levels, finest first, each 4x coarser than the previous one.""" 382 382 + g = global_level(GLOBAL_WIDTH).res 383 383 + resolutions = [max(proj.q0, native_km)] 384 384 + while resolutions[-1] * 4 < g / 1.5: 385 385 + resolutions.append(resolutions[-1] * 4) 386 386 + if resolutions[0] >= g / 1.5: 387 387 + return [] 388 388 + levels = [] 389 389 + for i, res in enumerate(resolutions): 390 390 + nxt = resolutions[i + 1] if i + 1 < len(resolutions) else g 391 391 + radius = proj.distance_for_demand(nxt) * 1.1 392 392 + if radius > 0: 393 393 + levels.append(local_level(proj, res, radius)) 394 394 + return levels 395 395 + 396 396 + 397 397 + def mipmaps(level: Level, q_max: float) -> list[Level]: 398 398 + out = [] 399 399 + cur = level 400 400 + while cur.res < q_max and cur.height >= 64 and cur.data is not None: 401 401 + d = cur.data.astype(np.uint16) 402 402 + d = (d[0::2, 0::2] + d[1::2, 0::2] + d[0::2, 1::2] + d[1::2, 1::2] + 2) // 4 403 403 + cur = Level( 404 404 + res=cur.res * 2, 405 405 + south=-90.0, 406 406 + west=-180.0, 407 407 + north=90.0, 408 408 + east=180.0, 409 409 + width=cur.width // 2, 410 410 + height=cur.height // 2, 411 411 + wrap=True, 412 412 + data=d.astype(np.uint8), 413 413 + ) 414 414 + out.append(cur) 415 415 + return out 416 416 + 417 417 + 418 418 + # -------------------------------------------------------------------------- 419 419 + # raster preparation: RGBA -> channels, last channel is data coverage, value 420 420 + # channels are premultiplied by coverage. Bump PREP_VERSION whenever a prep 421 421 + # function changes so cached prepared rasters are rebuilt from the raw ones. 422 422 + 423 423 + PREP_VERSION = 3 424 424 + 425 425 + 426 426 + class Colormap: 427 427 + """GIBS colormap: maps the rendered colours back to a 0..1 data value.""" 428 428 + 429 429 + def __init__(self, path: Path) -> None: 430 430 + import xml.etree.ElementTree as ET 431 431 + 432 432 + entries = [ 433 433 + e 434 434 + for e in ET.parse(path).iter() 435 435 + if e.tag.endswith("ColorMapEntry") 436 436 + and e.get("transparent") != "true" 437 437 + and e.get("nodata") != "true" 438 438 + ] 439 439 + if not entries: 440 440 + raise ValueError(f"{path}: no colormap entries") 441 441 + self.rgb = np.array([[int(c) for c in e.get("rgb", "").split(",")] for e in entries]) 442 442 + self.value = np.linspace(0.0, 1.0, len(entries)) 443 443 + 444 444 + def decode(self, rgba: np.ndarray) -> np.ndarray: 445 445 + out = np.zeros(rgba.shape[:2], np.float32) 446 446 + mask = rgba[..., 3] > 0 447 447 + px = rgba[mask][:, :3].astype(np.int64) 448 448 + if len(px) == 0: 449 449 + return out 450 450 + uniq, inv = np.unique((px[:, 0] << 16) | (px[:, 1] << 8) | px[:, 2], return_inverse=True) 451 451 + colors = np.stack([uniq >> 16, (uniq >> 8) & 255, uniq & 255], axis=-1) 452 452 + values = np.empty(len(uniq)) 453 453 + for i in range(0, len(uniq), 4096): 454 454 + d = ((colors[i : i + 4096, None, :] - self.rgb[None]) ** 2).sum(-1) 455 455 + values[i : i + 4096] = self.value[d.argmin(1)] 456 456 + out[mask] = values[inv] 457 457 + return out 458 458 + 459 459 + 460 460 + def to_u8(x: np.ndarray) -> np.ndarray: 461 461 + return np.round(np.clip(x, 0, 1) * 255).astype(np.uint8) 462 462 + 463 463 + 464 464 + def prep_rgb(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray: 465 465 + a = rgba[..., 3:4].astype(np.uint16) 466 466 + rgb = (rgba[..., :3].astype(np.uint16) * a + 127) // 255 467 467 + return np.concatenate([rgb.astype(np.uint8), rgba[..., 3:4]], axis=-1) 468 468 + 469 469 + 470 470 + def prep_lights(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray: 471 471 + """Black Marble -> city light intensity; drops its dim blue land texture.""" 472 472 + red = rgba[..., 0].astype(np.float32) / 255 473 473 + full = np.full(rgba.shape[:2], 255, np.uint8) 474 474 + return np.stack([to_u8(np.clip((red - 0.12) / 0.75, 0, 1) ** 1.2), full], axis=-1) 475 475 + 476 476 + 477 477 + def prep_cloud(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray: 478 478 + """Infrared 10.8 µm brightness (cold = bright) -> cloud opacity.""" 479 479 + v = rgba[..., 0].astype(np.float32) / 255 480 480 + cov = rgba[..., 3].astype(np.float32) / 255 481 481 + assert layer.transform is not None 482 482 + return np.stack([to_u8(layer.transform(v) * cov), rgba[..., 3]], axis=-1) 483 483 + 484 484 + 485 485 + def prep_scalar(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray: 486 486 + """Colormapped data value (or plain presence) -> opacity, full coverage.""" 487 487 + v = cmap.decode(rgba) if cmap else rgba[..., 3].astype(np.float32) / 255 488 488 + if layer.transform is not None: 489 489 + v = layer.transform(v) * (rgba[..., 3] > 0) 490 490 + full = np.full(rgba.shape[:2], 255, np.uint8) 491 491 + return np.stack([to_u8(v), full], axis=-1) 492 492 + 493 493 + 494 494 + def prep_tint(layer: "Layer", rgba: np.ndarray, cmap: Colormap | None) -> np.ndarray: 495 495 + """Colormapped value plus where it is defined, full coverage.""" 496 496 + assert cmap is not None 497 497 + present = (rgba[..., 3] > 0).astype(np.float32) 498 498 + full = np.full(rgba.shape[:2], 255, np.uint8) 499 499 + return np.stack([to_u8(cmap.decode(rgba) * present), to_u8(present), full], axis=-1) 500 500 + 501 501 + 502 502 + def combine_max(layers: list[np.ndarray]) -> np.ndarray: 503 503 + out = layers[0].copy() 504 504 + for older in layers[1:]: 505 505 + np.maximum(out, older, out=out) 506 506 + return out 507 507 + 508 508 + 509 509 + def combine_fill(layers: list[np.ndarray]) -> np.ndarray: 510 510 + """Fill where the newest has no data (channel -2 = presence) from older.""" 511 511 + out = layers[0].copy() 512 512 + for older in layers[1:]: 513 513 + gap = out[..., -2] == 0 514 514 + out[gap] = older[gap] 515 515 + return out 516 516 + 517 517 + 518 518 + # -------------------------------------------------------------------------- 519 519 + # fetching 520 520 + 521 521 + 522 522 + class FetchError(Exception): 523 523 + pass 524 524 + 525 525 + 526 526 + def http_get(url: str, timeout: float = 120) -> tuple[bytes, str]: 527 527 + last: Exception | None = None 528 528 + for attempt in range(3): 529 529 + try: 530 530 + req = urllib.request.Request(url, headers={"User-Agent": USER_AGENT}) 531 531 + with urllib.request.urlopen(req, timeout=timeout) as resp: 532 532 + return resp.read(), resp.headers.get("Content-Type", "") 533 533 + except Exception as e: # noqa: BLE001 - network errors come in many types 534 534 + last = e 535 535 + time.sleep(2 * (attempt + 1)) 536 536 + raise FetchError(f"{url}: {last}") 537 537 + 538 538 + 539 539 + def fetch_image(url: str) -> np.ndarray: 540 540 + body, ctype = http_get(url) 541 541 + if not ctype.startswith("image/"): 542 542 + raise FetchError(f"{url}: got {ctype}: {body[:200]!r}") 543 543 + return np.asarray(Image.open(io.BytesIO(body)).convert("RGBA")) 544 544 + 545 545 + 546 546 + def split_columns(level: Level) -> list[tuple[int, int]]: 547 547 + px = (level.east - level.west) / level.width 548 548 + cuts = {0, level.width} 549 549 + for meridian in (-180.0, 180.0): 550 550 + if level.west < meridian < level.east and not level.wrap: 551 551 + cuts.add(round((meridian - level.west) / px)) 552 552 + c = 0 553 553 + while c < level.width: 554 554 + cuts.add(c) 555 555 + c += MAX_TILE 556 556 + edges = sorted(x for x in cuts if 0 <= x <= level.width) 557 557 + return [(a, b) for a, b in zip(edges, edges[1:]) if b > a] 558 558 + 559 559 + 560 560 + def fetch_level(src: Wms, level: Level, t: str | None) -> np.ndarray: 561 561 + factor = max(1, min(src.oversample, 2 * MAX_LEVEL_SIDE // max(level.width, level.height))) 562 562 + big = dataclasses.replace(level, width=level.width * factor, height=level.height * factor, 563 563 + data=None) 564 564 + img = fetch_tiles(src, big, t) 565 565 + if factor > 1: 566 566 + img = np.asarray(Image.fromarray(img, "RGBA").reduce(factor)) 567 567 + return img 568 568 + 569 569 + 570 570 + def apply_gain(rgba: np.ndarray, gain: tuple[float, float, float] | None) -> np.ndarray: 571 571 + if gain is None: 572 572 + return rgba 573 573 + out = rgba.copy() 574 574 + out[..., :3] = to_u8(rgba[..., :3] / 255.0 * np.array(gain)) 575 575 + return out 576 576 + 577 577 + 578 578 + def fetch_tiles(src: Wms, level: Level, t: str | None) -> np.ndarray: 579 579 + px_lon = (level.east - level.west) / level.width 580 580 + px_lat = (level.north - level.south) / level.height 581 581 + rows = [(r, min(r + MAX_TILE, level.height)) for r in range(0, level.height, MAX_TILE)] 582 582 + img = np.zeros((level.height, level.width, 4), np.uint8) 583 583 + 584 584 + def one(r0: int, r1: int, c0: int, c1: int) -> None: 585 585 + w = level.west + c0 * px_lon 586 586 + e = level.west + c1 * px_lon 587 587 + mid = (w + e) / 2 588 588 + shift = -360.0 if mid > 180 else 360.0 if mid < -180 else 0.0 589 589 + w, e = max(-180.0, w + shift), min(180.0, e + shift) 590 590 + n = level.north - r0 * px_lat 591 591 + s = level.north - r1 * px_lat 592 592 + tile = fetch_image(src.url((s, w, n, e), c1 - c0, r1 - r0, t)) 593 593 + if tile.shape[:2] != (r1 - r0, c1 - c0): 594 594 + tile = np.asarray(Image.fromarray(tile).resize((c1 - c0, r1 - r0), Image.BILINEAR)) 595 595 + img[r0:r1, c0:c1] = tile 596 596 + 597 597 + with concurrent.futures.ThreadPoolExecutor(max_workers=TILE_WORKERS) as pool: 598 598 + futures = [pool.submit(one, r0, r1, c0, c1) 599 599 + for r0, r1 in rows for c0, c1 in split_columns(level)] 600 600 + for fut in futures: 601 601 + fut.result() # re-raise the first failure 602 602 + return img 603 603 + 604 604 + 605 605 + @dataclass 606 606 + class Layer: 607 607 + name: str 608 608 + local: Wms | None 609 609 + global_: Wms 610 610 + prep: Callable[["Layer", np.ndarray, Colormap | None], np.ndarray] 611 611 + combine: Callable[[list[np.ndarray]], np.ndarray] = combine_max 612 612 + colormap: str | None = None # GIBS colormap name 613 613 + transform: Callable[[np.ndarray], np.ndarray] | None = None 614 614 + required: bool = False 615 615 + cmap: Colormap | None = None 616 616 + levels: list[Level] = field(default_factory=list) 617 617 + log_res: np.ndarray | None = None 618 618 + 619 619 + 620 620 + class Cache: 621 621 + """Raw WMS rasters as PNG (refetched after the source TTL) and prepared 622 622 + rasters as .npy (rebuilt when a raw one changed or PREP_VERSION moved).""" 623 623 + 624 624 + def __init__(self, root: Path) -> None: 625 625 + self.root = root 626 626 + self.root.mkdir(parents=True, exist_ok=True) 627 627 + self.used: set[Path] = set() 628 628 + 629 629 + def _path(self, prefix: str, ident: object, suffix: str) -> Path: 630 630 + digest = hashlib.sha1(json.dumps(ident, sort_keys=True).encode()).hexdigest()[:16] 631 631 + path = self.root / f"{prefix}-{digest}{suffix}" 632 632 + self.used.add(path) 633 633 + return path 634 634 + 635 635 + def raw(self, src: Wms, level: Level, t: str | None) -> Path | None: 636 636 + ident = [src.layer, t, level.bbox, level.width, level.height, src.oversample] 637 637 + path = self._path(src.key, ident, ".png") 638 638 + if path.exists() and time.time() - path.stat().st_mtime < src.ttl: 639 639 + return path 640 640 + try: 641 641 + t0 = time.monotonic() 642 642 + img = fetch_level(src, level, t) 643 643 + tmp = path.with_name(f".{path.name}.tmp") 644 644 + Image.fromarray(img, "RGBA").save(tmp, format="PNG", compress_level=1) 645 645 + os.replace(tmp, path) 646 646 + log(f"fetched {src.key} {t or ''} {level.width}x{level.height} " 647 647 + f"@{level.res:.2f}km in {time.monotonic() - t0:.1f}s") 648 648 + except Exception as e: # noqa: BLE001 - any failure falls back to the cache 649 649 + log(f"warning: {src.key} {t or ''}: {e}" 650 650 + + ("; using cached copy" if path.exists() else "")) 651 651 + return path if path.exists() else None 652 652 + 653 653 + def get(self, layer: Layer, src: Wms, level: Level, when: dt.datetime) -> np.ndarray | None: 654 654 + raws = [p for t in src.times(when) if (p := self.raw(src, level, t)) is not None] 655 655 + if not raws: 656 656 + return None 657 657 + path = self._path( 658 658 + f"{layer.name}-prep", 659 659 + [PREP_VERSION, src.key, src.gain, [p.name for p in raws]], 660 660 + ".npy", 661 661 + ) 662 662 + newest = max(p.stat().st_mtime for p in raws) 663 663 + if not path.exists() or path.stat().st_mtime < newest: 664 664 + try: 665 665 + preps = [ 666 666 + layer.prep(layer, apply_gain(np.asarray(Image.open(p).convert("RGBA")), src.gain), 667 667 + layer.cmap) 668 668 + for p in raws 669 669 + ] 670 670 + tmp = path.with_name(f".{path.name}.tmp.npy") 671 671 + np.save(tmp, layer.combine(preps)) 672 672 + os.replace(tmp, path) 673 673 + except Exception as e: # noqa: BLE001 - a broken raster drops the level 674 674 + log(f"warning: {layer.name}: preparing {src.key}: {e}") 675 675 + return None 676 676 + return np.load(path, mmap_mode="r") 677 677 + 678 678 + def get_url(self, name: str, url: str, ttl: float) -> Path | None: 679 679 + path = self.root / name 680 680 + self.used.add(path) 681 681 + if path.exists() and time.time() - path.stat().st_mtime < ttl: 682 682 + return path 683 683 + try: 684 684 + body, _ = http_get(url, timeout=30) 685 685 + if not body: 686 686 + raise FetchError(f"{url}: empty response") 687 687 + tmp = path.with_name(f".{path.name}.tmp") 688 688 + tmp.write_bytes(body) 689 689 + os.replace(tmp, path) 690 690 + except Exception as e: # noqa: BLE001 691 691 + log(f"warning: {name}: {e}") 692 692 + return path if path.exists() else None 693 693 + 694 694 + def cleanup(self, max_age: float = 3 * 86400) -> None: 695 695 + for p in self.root.iterdir(): 696 696 + if p not in self.used and time.time() - p.stat().st_mtime > max_age: 697 697 + p.unlink() 698 698 + 699 699 + 700 700 + def snow_fraction(ndsi: np.ndarray) -> np.ndarray: 701 701 + # MODIS/VIIRS fractional snow cover from NDSI, ignoring the faint 702 702 + # detections that are mostly cloud edges and bright soil 703 703 + return smoothstep(0.35, 0.6, ndsi) 704 704 + 705 705 + 706 706 + def cloud_opacity(ir: np.ndarray) -> np.ndarray: 707 707 + return smoothstep(0.42, 0.78, ir) ** 1.2 708 708 + 709 709 + 710 710 + def layers() -> list[Layer]: 711 711 + day = 86400.0 712 712 + return [ 713 713 + Layer( 714 714 + "day", 715 715 + # s2cloudless is bluer than Blue Marble over the same land 716 716 + Wms("s2cloudless", EOX, "s2cloudless-2025", 0.01, 180 * day, "image/jpeg", 717 717 + oversample=4, gain=(1.0, 0.96, 0.85)), 718 718 + Wms("bluemarble", GIBS, "BlueMarble_NextGeneration", 0.5, 180 * day, "image/jpeg"), 719 719 + prep_rgb, 720 720 + required=True, 721 721 + ), 722 722 + Layer( 723 723 + "lights", 724 724 + None, 725 725 + Wms("blackmarble", GIBS, "VIIRS_Black_Marble", 0.5, 180 * day, "image/jpeg", 726 726 + time="2016-01-01"), 727 727 + prep_lights, 728 728 + ), 729 729 + Layer( 730 730 + "clouds", 731 731 + Wms("mtg-ir105", EUMETSAT, "mtg_fd:ir105_hrfi", 1.0, 600), 732 732 + Wms("georing-ir108", EUMETSAT, "mumi:worldcloudmap_ir108", 3.0, 3600), 733 733 + prep_cloud, 734 734 + transform=cloud_opacity, 735 735 + ), 736 736 + Layer( 737 737 + "snow", 738 738 + None, 739 739 + Wms("snow", GIBS, "VIIRS_NOAA20_NDSI_Snow_Cover", 0.375, 6 * 3600, 740 740 + days_back=(1, 2, 3)), 741 741 + prep_scalar, 742 742 + colormap="MODIS_NDSI_Snow_Cover", 743 743 + transform=snow_fraction, 744 744 + ), 745 745 + Layer( 746 746 + "ice", 747 747 + None, 748 748 + Wms("seaice", GIBS, "GHRSST_L4_MUR_Sea_Ice_Concentration", 1.0, 6 * 3600, 749 749 + days_back=(1, 2, 3)), 750 750 + prep_scalar, 751 751 + colormap="GHRSST_Sea_Ice_Concentration", 752 752 + transform=lambda c: smoothstep(0.15, 0.9, c), 753 753 + ), 754 754 + Layer( 755 755 + "fire", 756 756 + None, 757 757 + Wms("fire", GIBS, "VIIRS_NOAA20_Thermal_Anomalies_375m_All", 0.375, 3600, 758 758 + days_back=(0, 1)), 759 759 + prep_scalar, 760 760 + ), 761 761 + Layer( 762 762 + "chlorophyll", 763 763 + None, 764 764 + Wms("chlorophyll", GIBS, "VIIRS_NOAA20_Chlorophyll_a", 1.0, 6 * 3600, 765 765 + days_back=(1, 2, 3, 4)), 766 766 + prep_tint, 767 767 + combine=combine_fill, 768 768 + colormap="VIIRS_Chlorophyll", 769 769 + ), 770 770 + ] 771 771 + 772 772 + 773 773 + def load_layers(proj: LensMap, cache: Cache, when: dt.datetime) -> dict[str, Layer]: 774 774 + all_layers = [] 775 775 + for layer in layers(): 776 776 + if layer.colormap: 777 777 + path = cache.get_url(f"colormap-{layer.colormap}.xml", 778 778 + f"https://gibs.earthdata.nasa.gov/colormaps/v1.3/{layer.colormap}.xml", 779 779 + 30 * 86400) 780 780 + try: 781 781 + layer.cmap = Colormap(path) if path else None 782 782 + except Exception as e: # noqa: BLE001 783 783 + log(f"warning: colormap {layer.colormap}: {e}") 784 784 + if layer.cmap is None: 785 785 + log(f"warning: no colormap for {layer.name}, leaving it out") 786 786 + continue 787 787 + all_layers.append(layer) 788 788 + jobs = [] 789 789 + for layer in all_layers: 790 790 + local_src = layer.local or layer.global_ 791 791 + for lvl in plan_local(proj, local_src.native_km): 792 792 + jobs.append((layer, local_src, lvl)) 793 793 + jobs.append((layer, layer.global_, global_level(GLOBAL_WIDTH))) 794 794 + with concurrent.futures.ThreadPoolExecutor(max_workers=6) as pool: 795 795 + futures = [(layer, lvl, pool.submit(cache.get, layer, src, lvl, when)) 796 796 + for layer, src, lvl in jobs] 797 797 + for layer, lvl, fut in futures: 798 798 + lvl.data = fut.result() 799 799 + if lvl.data is not None: 800 800 + layer.levels.append(lvl) 801 801 + out = {} 802 802 + for layer in all_layers: 803 803 + if not layer.levels or not layer.levels[-1].wrap: 804 804 + if layer.required: 805 805 + raise SystemExit(f"no global data for required layer {layer.name}") 806 806 + log(f"warning: layer {layer.name} has no global data, leaving it out") 807 807 + continue 808 808 + layer.levels.extend(mipmaps(layer.levels[-1], proj.q_max)) 809 809 + # np.interp needs increasing resolutions; a size-capped local level 810 810 + # can end up coarser than the next one 811 811 + layer.log_res = np.maximum.accumulate(np.log([lvl.res for lvl in layer.levels])) 812 812 + out[layer.name] = layer 813 813 + return out 814 814 + 815 815 + 816 816 + # -------------------------------------------------------------------------- 817 817 + # sampling 818 818 + 819 819 + 820 820 + def sample_level(level: Level, lat: np.ndarray, lon: np.ndarray): 821 821 + data = level.data 822 822 + assert data is not None 823 823 + h, w = data.shape[:2] 824 824 + col = ((lon - level.west) % 360.0) / (level.east - level.west) * w - 0.5 825 825 + row = (level.north - lat) / (level.north - level.south) * h - 0.5 826 826 + if level.wrap: 827 827 + valid = np.ones(lat.shape, bool) 828 828 + c = np.floor(col) 829 829 + fc = (col - c)[:, None] 830 830 + c0 = c.astype(np.int64) % w 831 831 + c1 = (c0 + 1) % w 832 832 + else: 833 833 + valid = (col >= -0.5) & (col <= w - 0.5) & (row >= -0.5) & (row <= h - 0.5) 834 834 + col = np.clip(col, 0, w - 1) 835 835 + c0 = np.floor(col).astype(np.int64) 836 836 + c1 = np.minimum(c0 + 1, w - 1) 837 837 + fc = (col - c0)[:, None] 838 838 + row = np.clip(row, 0, h - 1) 839 839 + r0 = np.floor(row).astype(np.int64) 840 840 + r1 = np.minimum(r0 + 1, h - 1) 841 841 + fr = (row - r0)[:, None] 842 842 + top = data[r0, c0] * (1 - fc) + data[r0, c1] * fc 843 843 + bottom = data[r1, c0] * (1 - fc) + data[r1, c1] * fc 844 844 + return (top * (1 - fr) + bottom * fr) / 255.0, valid 845 845 + 846 846 + 847 847 + def sample(layer: Layer, lat: np.ndarray, lon: np.ndarray, logq: np.ndarray): 848 848 + """Trilinear sample: blend the two levels whose resolution brackets q. 849 849 + 850 850 + Returns premultiplied values and the total data coverage. 851 851 + """ 852 852 + n = len(layer.levels) 853 853 + t = np.interp(logq, layer.log_res, np.arange(n, dtype=np.float64)) 854 854 + channels = layer.levels[0].data.shape[-1] - 1 855 855 + out = np.zeros((len(lat), channels)) 856 856 + rem = np.ones(len(lat)) 857 857 + for i, lvl in enumerate(layer.levels): 858 858 + w = np.clip(i + 1 - t, 0.0, 1.0) 859 859 + idx = np.nonzero((w > 0) & (rem > 1e-4))[0] 860 860 + if len(idx) == 0: 861 861 + continue 862 862 + v, valid = sample_level(lvl, lat[idx], lon[idx]) 863 863 + wi = rem[idx] * w[idx] * valid 864 864 + out[idx] += v[:, :-1] * wi[:, None] 865 865 + rem[idx] -= wi * v[:, -1] 866 866 + return out, 1.0 - rem 867 867 + 868 868 + 869 869 + def lerp(a: np.ndarray, b, t: np.ndarray) -> np.ndarray: 870 870 + return a + (np.asarray(b) - a) * t[:, None] 871 871 + 872 872 + 873 873 + def subsolar_cos(lat, lon, sun_lat: float, sun_lon: float) -> np.ndarray: 874 874 + la, lo = np.radians(lat), np.radians(lon) 875 875 + sl, so = math.radians(sun_lat), math.radians(sun_lon) 876 876 + return np.sin(la) * math.sin(sl) + np.cos(la) * math.cos(sl) * np.cos(lo - so) 877 877 + 878 878 + 879 879 + # -------------------------------------------------------------------------- 880 880 + # cloud post-processing: the best live cloud data is ~1 km, far coarser than 881 881 + # the zoomed-in pixels around home. Synthesise the missing detail (noise 882 882 + # anchored to the ground, so it moves with the clouds, not the screen), and 883 883 + # light the clouds from the sun: relief shading plus shadows on the ground. 884 884 + 885 885 + CLOUD_DATA_KM = 1.0 # resolution of the finest cloud source 886 886 + CLOUD_HEIGHT_KM = 3.0 # assumed cloud height for the shadow offset 887 887 + 888 888 + 889 889 + def value_noise(p: np.ndarray) -> np.ndarray: 890 890 + """Smooth 3D value noise in [-1, 1] at integer lattice spacing.""" 891 891 + i = np.floor(p).astype(np.int64) 892 892 + f = p - i 893 893 + u = f * f * (3 - 2 * f) 894 894 + out = np.zeros(len(p)) 895 895 + for dx in (0, 1): 896 896 + for dy in (0, 1): 897 897 + for dz in (0, 1): 898 898 + h = ((i[:, 0] + dx) * 73856093) ^ ((i[:, 1] + dy) * 19349663) ^ ((i[:, 2] + dz) * 83492791) 899 899 + h = (h * 0x2545F491) & 0xFFFFFF 900 900 + w = (u[:, 0] if dx else 1 - u[:, 0]) * (u[:, 1] if dy else 1 - u[:, 1]) \ 901 901 + * (u[:, 2] if dz else 1 - u[:, 2]) 902 902 + out += w * (h / 0x7FFFFF - 1.0) 903 903 + return out 904 904 + 905 905 + 906 906 + def cloud_detail(a: np.ndarray, lat: np.ndarray, lon: np.ndarray, q: np.ndarray) -> np.ndarray: 907 907 + """Turn the soft, upsampled cloud edges into crisp fractal ones. 908 908 + 909 909 + Octaves (in km on the ground) are only added where they span at least 910 910 + two output px and lie below the data resolution, so the far, small parts 911 911 + of the map stay as they are. 912 912 + """ 913 913 + out = a.copy() 914 914 + flat_a, flat_q = a.ravel(), q.ravel() 915 915 + todo = np.nonzero((flat_a > 0.02) & (flat_a < 0.98) & (flat_q < 2 * CLOUD_DATA_KM))[0] 916 916 + octaves = [(0.4, 0.35), (0.8, 0.5), (1.6, 0.75), (3.2, 1.0)] # wavelength km, amplitude 917 917 + total = sum(amp for _, amp in octaves) 918 918 + flat_out = out.ravel() 919 919 + for c0 in range(0, len(todo), 1_000_000): 920 920 + idx = todo[c0 : c0 + 1_000_000] 921 921 + la, lo = np.radians(lat.ravel()[idx]), np.radians(lon.ravel()[idx]) 922 922 + p = R_EARTH * np.stack([np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)], -1) 923 923 + qi, ai = flat_q[idx], flat_a[idx] 924 924 + n = np.zeros(len(idx)) 925 925 + weight = np.zeros(len(idx)) 926 926 + for wavelength, amp in octaves: 927 927 + w = smoothstep(2.0, 4.0, wavelength / qi) * amp 928 928 + n += w * value_noise(p / wavelength) 929 929 + weight += w 930 930 + detail = np.clip(weight / total * 1.5, 0, 1) * smoothstep(2 * CLOUD_DATA_KM, 0.5 * CLOUD_DATA_KM, qi) 931 931 + nn = n / np.maximum(weight, 1e-6) # fBm in [-1, 1] 932 932 + # modulate rather than threshold: thin parts get the most texture, 933 933 + # dense cores only some, and translucency survives 934 934 + textured = np.clip(ai * (1 + 0.9 * nn * (1 - 0.6 * ai)), 0, 1) 935 935 + textured = textured + (smoothstep(0.15, 0.85, textured) - textured) * 0.5 936 936 + flat_out[idx] = ai + (textured - ai) * detail 937 937 + return out 938 938 + 939 939 + 940 940 + def box_blur(a: np.ndarray, r: int) -> np.ndarray: 941 941 + """Approximate Gaussian blur (three box passes per axis).""" 942 942 + if r < 1: 943 943 + return a 944 944 + out = a.astype(np.float32) 945 945 + for axis in (0, 1): 946 946 + for _ in range(3): 947 947 + pad = [(0, 0), (0, 0)] 948 948 + pad[axis] = (r + 1, r) 949 949 + c = np.cumsum(np.pad(out, pad, mode="edge"), axis=axis, dtype=np.float64) 950 950 + hi = np.take(c, np.arange(2 * r + 1, c.shape[axis]), axis=axis) 951 951 + lo = np.take(c, np.arange(0, c.shape[axis] - 2 * r - 1), axis=axis) 952 952 + out = ((hi - lo) / (2 * r + 1)).astype(np.float32) 953 953 + return out 954 954 + 955 955 + 956 956 + def gather(a: np.ndarray, x: np.ndarray, y: np.ndarray) -> np.ndarray: 957 957 + """Bilinear lookup of a 2D array at fractional pixel positions.""" 958 958 + h, w = a.shape 959 959 + x = np.clip(x, 0, w - 1.001) 960 960 + y = np.clip(y, 0, h - 1.001) 961 961 + x0, y0 = np.floor(x).astype(np.int64), np.floor(y).astype(np.int64) 962 962 + fx, fy = x - x0, y - y0 963 963 + top = a[y0, x0] * (1 - fx) + a[y0, x0 + 1] * fx 964 964 + bottom = a[y0 + 1, x0] * (1 - fx) + a[y0 + 1, x0 + 1] * fx 965 965 + return top * (1 - fy) + bottom * fy 966 966 + 967 967 + 968 968 + def render_base(proj: LensMap, lay: dict[str, Layer], sun: tuple[float, float]) -> np.ndarray: 969 969 + W, H = proj.width, proj.height 970 970 + ground = np.zeros((H, W, 3), np.float32) 971 971 + cloud = np.zeros((H, W), np.float32) 972 972 + cosz_f = np.zeros((H, W), np.float32) 973 973 + lights = np.zeros((H, W), np.float32) 974 974 + lat_f = np.zeros((H, W), np.float32) 975 975 + lon_f = np.zeros((H, W), np.float32) 976 976 + q_f = np.zeros((H, W), np.float32) 977 977 + chunk = max(1, 1_000_000 // W) 978 978 + xs = np.arange(W) + 0.5 979 979 + snow_rgb = np.array([0.87, 0.87, 0.87]) 980 980 + ice_rgb = np.array([0.90, 0.94, 1.00]) 981 981 + fire_rgb = np.array([0.94, 0.52, 0.05]) 982 982 + # chlorophyll: turquoise for moderately productive water, green for blooms 983 983 + plankton_rgb = np.array([0.06, 0.32, 0.40]) 984 984 + bloom_rgb = np.array([0.14, 0.40, 0.16]) 985 985 + # the night side is the (high resolution) day map in faint moonlight 986 986 + moonlight = np.array([0.20, 0.24, 0.36]) 987 987 + for y0 in range(0, H, chunk): 988 988 + ys = np.arange(y0, min(H, y0 + chunk)) + 0.5 989 989 + rows = slice(y0, y0 + len(ys)) 990 990 + X, Y = np.meshgrid(xs, ys) 991 991 + lat, lon, q = proj.inverse(X.ravel(), Y.ravel()) 992 992 + logq = np.log(q) 993 993 + 994 994 + def get(name: str): 995 995 + return sample(lay[name], lat, lon, logq) if name in lay else None 996 996 + 997 997 + day, cov = get("day") 998 998 + day /= np.maximum(cov, 1e-6)[:, None] 999 999 + 1000 1000 + if (s := get("chlorophyll")) is not None: 1001 1001 + # the value is log-scaled 0.01..20 mg/m³: 0.30 is ~0.1 mg/m³ 1002 1002 + # (productive shelf water), 0.75 is ~3 mg/m³ (a bloom) 1003 1003 + v, _ = s 1004 1004 + present = v[:, 1] 1005 1005 + value = v[:, 0] / np.maximum(present, 1e-6) 1006 1006 + tint = lerp(np.tile(plankton_rgb, (len(value), 1)), bloom_rgb, 1007 1007 + smoothstep(0.55, 0.85, value)) 1008 1008 + day = lerp(day, tint, present * smoothstep(0.30, 0.75, value) * 0.7) 1009 1009 + if (s := get("ice")) is not None: 1010 1010 + day = lerp(day, ice_rgb, s[0][:, 0] * 0.9) 1011 1011 + if (s := get("snow")) is not None: 1012 1012 + day = lerp(day, snow_rgb, s[0][:, 0] * 0.85) 1013 1013 + 1014 1014 + cosz = subsolar_cos(lat, lon, *sun) 1015 1015 + light = smoothstep(-0.10, 0.08, cosz) 1016 1016 + day_lit = day * (0.55 + 0.45 * np.sqrt(np.clip(cosz, 0, 1)))[:, None] 1017 1017 + color = lerp(day * moonlight, day_lit, light) 1018 1018 + if (s := get("fire")) is not None: 1019 1019 + color = np.maximum(color, fire_rgb * 0.8 * s[0][:, 0][:, None]) 1020 1020 + if (s := get("clouds")) is not None: 1021 1021 + cloud[rows] = np.clip(s[0][:, 0], 0, 1).reshape(len(ys), W) 1022 1022 + if (s := get("lights")) is not None: 1023 1023 + lights[rows] = (s[0][:, 0] * (1 - light)).reshape(len(ys), W) 1024 1024 + 1025 1025 + ground[rows] = color.reshape(len(ys), W, 3) 1026 1026 + cosz_f[rows] = cosz.reshape(len(ys), W) 1027 1027 + lat_f[rows] = lat.reshape(len(ys), W) 1028 1028 + lon_f[rows] = lon.reshape(len(ys), W) 1029 1029 + q_f[rows] = q.reshape(len(ys), W) 1030 1030 + 1031 1031 + scale = H / 1800 1032 1032 + light = smoothstep(-0.10, 0.08, cosz_f) 1033 1033 + if "clouds" in lay: 1034 1034 + cloud = cloud_detail(cloud, lat_f, lon_f, q_f) 1035 1035 + del lat_f, lon_f 1036 1036 + 1037 1037 + # sun direction on screen: cos(zenith) grows towards the sun 1038 1038 + gy, gx = np.gradient(cosz_f) 1039 1039 + norm = np.hypot(gx, gy) 1040 1040 + lx = np.where(norm > 1e-9, gx / np.maximum(norm, 1e-9), 0.0).astype(np.float32) 1041 1041 + ly = np.where(norm > 1e-9, gy / np.maximum(norm, 1e-9), 0.0).astype(np.float32) 1042 1042 + 1043 1043 + # shadows: the cloud seen from the ground towards the sun, displaced by 1044 1044 + # height / tan(elevation), only where the sun is up 1045 1045 + elev = np.arcsin(np.clip(cosz_f, 0.0, 1.0)) 1046 1046 + offset = np.clip(CLOUD_HEIGHT_KM / np.maximum(np.tan(elev), 0.1) / q_f, 0, 40 * scale) 1047 1047 + yy, xx = np.mgrid[0:H, 0:W].astype(np.float32) 1048 1048 + soft = box_blur(cloud, max(1, round(1.5 * scale))) 1049 1049 + shadow = gather(soft, xx + lx * offset, yy + ly * offset) 1050 1050 + del xx, yy 1051 1051 + ground *= (1 - 0.45 * shadow * light * smoothstep(0.0, 0.15, cosz_f))[..., None] 1052 1052 + 1053 1053 + # relief: treat the (blurred) opacity as cloud height and light its 1054 1054 + # slopes; a slope facing the sun has its height falling towards the sun 1055 1055 + height = box_blur(cloud, max(1, round(2 * scale))) 1056 1056 + hy, hx = np.gradient(height) 1057 1057 + relief = np.clip(1 - 6.0 * (hx * lx + hy * ly) / scale, 0.6, 1.12) 1058 1058 + 1059 1059 + sunlit = light * (0.6 + 0.4 * np.sqrt(np.clip(cosz_f, 0, 1))) * relief 1060 1060 + cloud_rgb = ( 1061 1061 + np.array([0.11, 0.12, 0.16], np.float32) * (1 - light)[..., None] 1062 1062 + + np.array([0.88, 0.89, 0.90], np.float32) * sunlit[..., None] 1063 1063 + ) 1064 1064 + out = ground + (cloud_rgb - ground) * cloud[..., None] 1065 1065 + 1066 1066 + # city lights as a glow drawn in screen space: smooth at any zoom, so the 1067 1067 + # ~1 km squares of the source never show; clouds dim them 1068 1068 + lights *= 1 - 0.75 * cloud 1069 1069 + img = Image.fromarray(to_u8(lights), "L") 1070 1070 + core = np.asarray(img.filter(ImageFilter.GaussianBlur(0.9 * scale)), np.float32) / 255 1071 1071 + halo = np.asarray(img.filter(ImageFilter.GaussianBlur(6 * scale)), np.float32) / 255 1072 1072 + glow = np.minimum(core * 1.1 + halo * 0.8, 1.5)[..., None] 1073 1073 + out += glow * np.array([1.0, 0.80, 0.50], np.float32) 1074 1074 + return to_u8(out) 1075 1075 + 1076 1076 + 1077 1077 + # -------------------------------------------------------------------------- 1078 1078 + # sky: sun, moon, planets and constellations at their sub-points 1079 1079 + 1080 1080 + 1081 1081 + PLANETS = { 1082 1082 + "mercury": "☿", 1083 1083 + "venus": "♀", 1084 1084 + "mars": "♂", 1085 1085 + "jupiter": "♃", 1086 1086 + "saturn": "♄", 1087 1087 + "uranus": "♅", 1088 1088 + "neptune": "♆", 1089 1089 + } 1090 1090 + 1091 1091 + 1092 1092 + @dataclass 1093 1093 + class Sky: 1094 1094 + bodies: dict[str, tuple[float, float]] 1095 1095 + moon_elongation: float 1096 1096 + lines: list[tuple[tuple[float, float], tuple[float, float]]] 1097 1097 + 1098 1098 + 1099 1099 + def wrap180(x): 1100 1100 + return (np.asarray(x) + 180.0) % 360.0 - 180.0 1101 1101 + 1102 1102 + 1103 1103 + def compute_sky(when: dt.datetime, constellations: Path | None) -> Sky: 1104 1104 + warnings.filterwarnings("ignore", module="astropy") 1105 1105 + warnings.filterwarnings("ignore", module="erfa") 1106 1106 + from astropy.utils import iers 1107 1107 + 1108 1108 + # offline: bundled IERS tables, UT1-UTC drift is far below a pixel 1109 1109 + iers.conf.auto_download = False 1110 1110 + iers.conf.auto_max_age = None 1111 1111 + iers.conf.iers_degraded_accuracy = "warn" 1112 1112 + import astropy.units as u 1113 1113 + from astropy.coordinates import TETE, SkyCoord, get_body 1114 1114 + from astropy.time import Time 1115 1115 + 1116 1116 + t = Time(when) 1117 1117 + gast = t.sidereal_time("apparent", "greenwich").deg 1118 1118 + frame = TETE(obstime=t) 1119 1119 + 1120 1120 + def subpoints(coord): 1121 1121 + c = coord.transform_to(frame) 1122 1122 + return np.asarray(c.dec.deg), wrap180(np.asarray(c.ra.deg) - gast) 1123 1123 + 1124 1124 + bodies = {} 1125 1125 + raw = {} 1126 1126 + for name in ["sun", "moon", *PLANETS]: 1127 1127 + raw[name] = get_body(name, t) 1128 1128 + la, lo = subpoints(raw[name]) 1129 1129 + bodies[name] = (float(la), float(lo)) 1130 1130 + elong = float(raw["sun"].separation(raw["moon"]).deg) 1131 1131 + 1132 1132 + lines = [] 1133 1133 + if constellations is not None: 1134 1134 + feats = json.loads(constellations.read_text())["features"] 1135 1135 + strips = [ 1136 1136 + line for f in feats for line in f["geometry"]["coordinates"] if len(line) > 1 1137 1137 + ] 1138 1138 + pts = np.array([p for line in strips for p in line], dtype=np.float64) 1139 1139 + la, lo = subpoints(SkyCoord(ra=pts[:, 0] * u.deg, dec=pts[:, 1] * u.deg, frame="icrs")) 1140 1140 + i = 0 1141 1141 + for line in strips: 1142 1142 + for k in range(len(line) - 1): 1143 1143 + lines.append(((la[i + k], lo[i + k]), (la[i + k + 1], lo[i + k + 1]))) 1144 1144 + i += len(line) 1145 1145 + return Sky(bodies, elong, lines) 1146 1146 + 1147 1147 + 1148 1148 + def to_vec(lat: float, lon: float) -> np.ndarray: 1149 1149 + la, lo = math.radians(lat), math.radians(lon) 1150 1150 + return np.array([math.cos(la) * math.cos(lo), math.cos(la) * math.sin(lo), math.sin(la)]) 1151 1151 + 1152 1152 + 1153 1153 + def great_circle(a, b, proj: LensMap, max_px: float = 3.0) -> list[list[tuple[float, float]]]: 1154 1154 + """Project the great circle arc a->b, split where it wraps around the map.""" 1155 1155 + va, vb = to_vec(*a), to_vec(*b) 1156 1156 + omega = math.acos(max(-1.0, min(1.0, float(va @ vb)))) 1157 1157 + # sample densely, then drop points that add nothing on screen 1158 1158 + n = 2 1159 1159 + while n < 4096: 1160 1160 + ts = np.linspace(0, 1, n) 1161 1161 + if omega < 1e-9: 1162 1162 + vs = np.outer(1 - ts, va) + np.outer(ts, vb) 1163 1163 + else: 1164 1164 + vs = (np.outer(np.sin((1 - ts) * omega), va) + np.outer(np.sin(ts * omega), vb)) / math.sin(omega) 1165 1165 + lat = np.degrees(np.arcsin(np.clip(vs[:, 2], -1, 1))) 1166 1166 + lon = np.degrees(np.arctan2(vs[:, 1], vs[:, 0])) 1167 1167 + x, y = proj.forward(lat, lon) 1168 1168 + step = np.hypot(np.diff(x), np.diff(y)) 1169 1169 + if step.max() <= max_px: 1170 1170 + break 1171 1171 + n *= 2 1172 1172 + pieces, cur = [], [(x[0], y[0])] 1173 1173 + for k, s in enumerate(step): 1174 1174 + if s > proj.width * 0.25: 1175 1175 + if len(cur) > 1: 1176 1176 + pieces.append(cur) 1177 1177 + cur = [] 1178 1178 + cur.append((x[k + 1], y[k + 1])) 1179 1179 + if len(cur) > 1: 1180 1180 + pieces.append(cur) 1181 1181 + return pieces 1182 1182 + 1183 1183 + 1184 1184 + def sun_icon(path: Path | None, size: int) -> Image.Image: 1185 1185 + if path is not None: 1186 1186 + try: 1187 1187 + img = Image.open(path).convert("RGB").resize((size, size), Image.LANCZOS) 1188 1188 + yy, xx = np.mgrid[:size, :size] 1189 1189 + r = np.hypot(xx - size / 2 + 0.5, yy - size / 2 + 0.5) / (size / 2) 1190 1190 + alpha = np.clip((0.86 - r) / 0.06, 0, 1) 1191 1191 + rgba = np.dstack([np.asarray(img), (alpha * 255).astype(np.uint8)]) 1192 1192 + return Image.fromarray(rgba, "RGBA") 1193 1193 + except Exception as e: # noqa: BLE001 1194 1194 + log(f"warning: sun image: {e}") 1195 1195 + img = Image.new("RGBA", (size, size)) 1196 1196 + ImageDraw.Draw(img).ellipse((size * 0.2, size * 0.2, size * 0.8, size * 0.8), fill=(255, 210, 80, 255)) 1197 1197 + return img 1198 1198 + 1199 1199 + 1200 1200 + def moon_icon(size: int, elongation: float, sun_dir: tuple[float, float]) -> Image.Image: 1201 1201 + """Moon disc lit from screen direction sun_dir with the real phase.""" 1202 1202 + phase = math.radians(180.0 - elongation) # sun-moon-earth angle 1203 1203 + ex, ey = sun_dir 1204 1204 + yy, xx = np.mgrid[:size, :size] 1205 1205 + x = (xx + 0.5 - size / 2) / (size / 2) 1206 1206 + y = -(yy + 0.5 - size / 2) / (size / 2) 1207 1207 + r2 = x * x + y * y 1208 1208 + z = np.sqrt(np.clip(1 - r2, 0, 1)) 1209 1209 + lit = (x * ex + y * ey) * math.sin(phase) + z * math.cos(phase) 1210 1210 + shade = np.clip(lit * 4, 0, 1) * 0.85 + 0.12 1211 1211 + alpha = np.clip((1 - np.sqrt(r2)) * size / 2, 0, 1) 1212 1212 + g = (shade * 235).astype(np.uint8) 1213 1213 + rgba = np.dstack([g, g, (shade * 225).astype(np.uint8), (alpha * 255).astype(np.uint8)]) 1214 1214 + return Image.fromarray(rgba, "RGBA") 1215 1215 + 1216 1216 + 1217 1217 + def paste_centered(canvas: Image.Image, icon: Image.Image, x: float, y: float) -> None: 1218 1218 + canvas.alpha_composite(icon, (int(round(x - icon.width / 2)), int(round(y - icon.height / 2)))) 1219 1219 + 1220 1220 + 1221 1221 + def draw_lines(canvas: Image.Image, proj: LensMap, sky: Sky, scale: float) -> None: 1222 1222 + layer = Image.new("RGBA", canvas.size) 1223 1223 + d = ImageDraw.Draw(layer) 1224 1224 + width = max(1, round(1.2 * scale)) 1225 1225 + for a, b in sky.lines: 1226 1226 + for piece in great_circle(a, b, proj): 1227 1227 + d.line(piece, fill=(255, 255, 255, 110), width=width, joint="curve") 1228 1228 + canvas.alpha_composite(layer) 1229 1229 + 1230 1230 + 1231 1231 + def draw_bodies(canvas, proj, sky: Sky, scale: float, font_path, sun_image) -> None: 1232 1232 + sx, sy = proj.forward(*sky.bodies["sun"]) 1233 1233 + mx, my = proj.forward(*sky.bodies["moon"]) 1234 1234 + font = ImageFont.truetype(str(font_path), max(8, round(34 * scale))) if font_path else None 1235 1235 + d = ImageDraw.Draw(canvas) 1236 1236 + for name, glyph in PLANETS.items(): 1237 1237 + x, y = proj.forward(*sky.bodies[name]) 1238 1238 + if font is not None: 1239 1239 + d.text((float(x), float(y)), glyph, font=font, anchor="mm", fill=(254, 128, 25, 255), 1240 1240 + stroke_width=max(1, round(2 * scale)), stroke_fill=(0, 0, 0, 200)) 1241 1241 + else: 1242 1242 + r = 6 * scale 1243 1243 + d.ellipse((x - r, y - r, x + r, y + r), fill=(254, 128, 25, 255)) 1244 1244 + dx, dy = float(sx - mx), float(my - sy) 1245 1245 + norm = math.hypot(dx, dy) or 1.0 1246 1246 + paste_centered(canvas, moon_icon(max(8, round(50 * scale)), sky.moon_elongation, (dx / norm, dy / norm)), 1247 1247 + float(mx), float(my)) 1248 1248 + paste_centered(canvas, sun_icon(sun_image, max(8, round(100 * scale))), float(sx), float(sy)) 1249 1249 + 1250 1250 + 1251 1251 + def draw_markers(canvas, proj, markers: list[tuple[float, float]], icon_path, scale: float) -> None: 1252 1252 + if not markers or icon_path is None: 1253 1253 + return 1254 1254 + size = max(8, round(20 * scale)) 1255 1255 + icon = Image.open(icon_path).convert("RGBA").resize((size, size), Image.LANCZOS) 1256 1256 + for lat, lon in markers: 1257 1257 + x, y = proj.forward(lat, lon) 1258 1258 + paste_centered(canvas, icon, float(x), float(y)) 1259 1259 + 1260 1260 + 1261 1261 + def load_markers(path: Path | None) -> list[tuple[float, float]]: 1262 1262 + if path is None: 1263 1263 + return [] 1264 1264 + try: 1265 1265 + data = json.loads(path.read_text()) 1266 1266 + except Exception as e: # noqa: BLE001 1267 1267 + log(f"warning: markers: {e}") 1268 1268 + return [] 1269 1269 + entries = data.values() if isinstance(data, dict) else data 1270 1270 + out = [] 1271 1271 + for v in entries: 1272 1272 + if isinstance(v, dict) and v.get("latitude") is not None and v.get("longitude") is not None: 1273 1273 + out.append((float(v["latitude"]), float(v["longitude"]))) 1274 1274 + return out 1275 1275 + 1276 1276 + 1277 1277 + def save(img: Image.Image, path: Path, **kw) -> None: 1278 1278 + tmp = path.with_name(f".{path.name}.tmp") 1279 1279 + img.save(tmp, format=path.suffix.lstrip(".").replace("jpg", "jpeg").upper(), **kw) 1280 1280 + os.replace(tmp, path) 1281 1281 + 1282 1282 + 1283 1283 + # -------------------------------------------------------------------------- 1284 1284 + 1285 1285 + 1286 1286 + def main() -> None: 1287 1287 + env = os.environ.get 1288 1288 + p = argparse.ArgumentParser(prog="realwallpaper", description=__doc__.splitlines()[0]) 1289 1289 + p.add_argument("--dir", type=Path, default=Path("."), 1290 1290 + help="output and cache directory (default: .)") 1291 1291 + p.add_argument("--lat", type=float, default=52.520008, 1292 1292 + help="home latitude in degrees, south negative (default Berlin)") 1293 1293 + p.add_argument("--lon", type=float, default=13.404954, 1294 1294 + help="home longitude in degrees, west negative (default Berlin)") 1295 1295 + p.add_argument("--width", type=int, default=3840) 1296 1296 + p.add_argument("--height", type=int, default=2160) 1297 1297 + p.add_argument("--magnification", type=float, default=6.0, 1298 1298 + help="zoom at home relative to the plain map (default 6)") 1299 1299 + p.add_argument("--lens", type=float, default=1.0, 1300 1300 + help="lens size relative to the screen; 1 reaches the edges (default 1)") 1301 1301 + p.add_argument("--lens-shape", choices=["circle", "ellipse", "rect"], default="ellipse", 1302 1302 + help="circle; ellipse or rect follow the screen's aspect (default ellipse)") 1303 1303 + p.add_argument("--focus", type=float, default=0.0, 1304 1304 + help="share of the lens zoomed undistorted before the transition (default 0)") 1305 1305 + p.add_argument("--falloff", type=float, default=0.65, 1306 1306 + help="where across the lens (0..1) the distortion peaks; higher keeps " 1307 1307 + "the middle closer to a plain zoom (default 0.65)") 1308 1308 + p.add_argument("--projection", choices=["lens", "plain", "mercator"], default="lens", 1309 1309 + help="lens: centred on home with the zoom lens; plain: centred on home " 1310 1310 + "without distortion; mercator: classic north-up map (default lens)") 1311 1311 + p.add_argument("--marker-url", default=None, 1312 1312 + help="JSON with {name: {latitude, longitude}} for krebs markers") 1313 1313 + p.add_argument("--time", type=dt.datetime.fromisoformat, default=None, 1314 1314 + help="render for this UTC time instead of now (ISO 8601)") 1315 1315 + p.add_argument("--no-archive", action="store_true", help="do not keep a jpg in archive/") 1316 1316 + args = p.parse_args() 1317 1317 + 1318 1318 + when = args.time or dt.datetime.now(dt.timezone.utc) 1319 1319 + if when.tzinfo is None: 1320 1320 + when = when.replace(tzinfo=dt.timezone.utc) 1321 1321 + when = when.astimezone(dt.timezone.utc) 1322 1322 + 1323 1323 + out_dir = args.dir.resolve() 1324 1324 + out_dir.mkdir(parents=True, exist_ok=True) 1325 1325 + cache = Cache(out_dir / "cache") 1326 1326 + proj = LensMap(args.lat, args.lon, args.width, args.height, 1327 1327 + args.magnification, args.lens, args.lens_shape, args.focus, args.falloff, 1328 1328 + args.projection) 1329 1329 + home_x, home_y = proj.forward(args.lat, args.lon) 1330 1330 + home_q = proj.inverse(np.array([home_x]), np.array([home_y]))[2][0] 1331 1331 + log(f"projection: {args.projection}, home at {home_q * 1000:.0f} m/px, " 1332 1332 + f"finest {proj.q0 * 1000:.0f} m/px") 1333 1333 + 1334 1334 + font = env("REALWALLPAPER_FONT") 1335 1335 + constellations = env("REALWALLPAPER_CONSTELLATIONS") 1336 1336 + krebs_icon = env("REALWALLPAPER_KREBS_ICON") 1337 1337 + 1338 1338 + t0 = time.monotonic() 1339 1339 + sun_image = cache.get_url("sun.jpg", "https://sdo.gsfc.nasa.gov/assets/img/latest/latest_512_0171.jpg", 3600) 1340 1340 + marker_file = cache.get_url("markers.json", args.marker_url, 0) if args.marker_url else None 1341 1341 + lay = load_layers(proj, cache, when) 1342 1342 + log(f"data ready in {time.monotonic() - t0:.1f}s") 1343 1343 + 1344 1344 + try: 1345 1345 + sky = compute_sky(when, Path(constellations) if constellations else None) 1346 1346 + except Exception as e: # noqa: BLE001 - the map is still useful without the sky 1347 1347 + log(f"warning: sky: {e}") 1348 1348 + sky = None 1349 1349 + sun = sky.bodies["sun"] if sky else (0.0, 0.0) 1350 1350 + 1351 1351 + t0 = time.monotonic() 1352 1352 + base = Image.fromarray(render_base(proj, lay, sun)).convert("RGBA") 1353 1353 + log(f"rendered in {time.monotonic() - t0:.1f}s") 1354 1354 + 1355 1355 + scale = args.height / 1800 1356 1356 + markers = load_markers(marker_file) 1357 1357 + 1358 1358 + plain = base 1359 1359 + marker = plain.copy() 1360 1360 + if sky: 1361 1361 + draw_bodies(marker, proj, sky, scale, font, sun_image) 1362 1362 + krebs = marker.copy() 1363 1363 + draw_markers(krebs, proj, markers, krebs_icon, scale) 1364 1364 + stars = plain.copy() 1365 1365 + if sky: 1366 1366 + draw_lines(stars, proj, sky, scale) 1367 1367 + draw_bodies(stars, proj, sky, scale, font, sun_image) 1368 1368 + draw_markers(stars, proj, markers, krebs_icon, scale) 1369 1369 + 1370 1370 + for name, img in [ 1371 1371 + ("realwallpaper.png", plain), 1372 1372 + ("realwallpaper-marker.png", marker), 1373 1373 + ("realwallpaper-krebs.png", krebs), 1374 1374 + ("realwallpaper-krebs-stars.png", stars), 1375 1375 + ]: 1376 1376 + save(img.convert("RGB"), out_dir / name, compress_level=3) 1377 1377 + if not args.no_archive: 1378 1378 + (out_dir / "archive").mkdir(exist_ok=True) 1379 1379 + save(krebs.convert("RGB"), out_dir / "archive" / f"{when.strftime('%Y-%m-%dT%H:%M:%SZ')}.jpg", 1380 1380 + quality=88) 1381 1381 + cache.cleanup() 1382 1382 + log("done") 1383 1383 + 1384 1384 + 1385 1385 + if __name__ == "__main__": 1386 1386 + main()

+43

tools/realwallpaper/usage.kdl

View file Reviewed

··· 1 1 + bin "realwallpaper" 2 2 + about "Render a live world map wallpaper: oblique Mercator centred on home with a fisheye lens (clouds, day/night, snow, ice, fire, sky)" 3 3 + 4 4 + flag "--dir" help="Output and cache directory (default .)" { 5 5 + arg "<dir>" 6 6 + } 7 7 + flag "--lat" help="Home latitude (default Berlin)" { 8 8 + arg "<deg>" 9 9 + } 10 10 + flag "--lon" help="Home longitude (default Berlin)" { 11 11 + arg "<deg>" 12 12 + } 13 13 + flag "--width" help="Output width in px (default 3840)" { 14 14 + arg "<px>" 15 15 + } 16 16 + flag "--height" help="Output height in px (default 2160)" { 17 17 + arg "<px>" 18 18 + } 19 19 + flag "--magnification" help="Zoom at home relative to the plain map (default 6)" { 20 20 + arg "<factor>" 21 21 + } 22 22 + flag "--lens" help="Lens size relative to the screen; 1 reaches the edges (default 1)" { 23 23 + arg "<fraction>" 24 24 + } 25 25 + flag "--lens-shape" help="circle; ellipse or rect follow the screen's aspect (default ellipse)" { 26 26 + arg "<shape>" 27 27 + } 28 28 + flag "--focus" help="Share of the lens zoomed undistorted before the transition (default 0)" { 29 29 + arg "<fraction>" 30 30 + } 31 31 + flag "--falloff" help="Where across the lens (0..1) the distortion peaks; higher keeps the middle closer to a plain zoom (default 0.65)" { 32 32 + arg "<position>" 33 33 + } 34 34 + flag "--projection" help="lens: centred on home with the zoom lens; plain: centred on home without distortion; mercator: classic north-up map (default lens)" { 35 35 + arg "<mode>" 36 36 + } 37 37 + flag "--marker-url" help="JSON {name: {latitude, longitude}} for krebs markers" { 38 38 + arg "<url>" 39 39 + } 40 40 + flag "--time" help="Render for this UTC time instead of now (ISO 8601)" { 41 41 + arg "<time>" 42 42 + } 43 43 + flag "--no-archive" help="Do not keep a jpg copy in archive/"