Geothermal Animated

This was derived from NCEP Reanalysis data, in the tradition of Measuring Geothermal …

Enjoy 🙂 -Zoe

Addendum

geochg.sh:

# source geochg.sh
# Zoe Phin 2020/03/13
    
F=(0 ulwrf dswrf uswrf lhtfl shtfl)                                                  
O=(0 3201.5 3086.5 3131.5 856.5 2176.5)

require() { sudo apt install nco gnuplot imagemagick; } # Linux Only
    
download() {
    b="ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis2.derived/gaussian_grid"
    for i in ${F[*]}; do wget -O $i.nc -c $b/$i.sfc.mon.mean.nc; done
}

extract() {
    for t in {000..491}; do echo "$t" >&2
        for i in {1..5}; do 
            ncks --trd -HC ${F[$i]}.nc -v ${F[$i]} -d time,$t | sed \$d | awk -F[=\ ] -vO=${O[$i]} '{ 
                printf "%7s %7s %7.3f\n", $4, $6, $8/10+O }' > .f$i
        done
        paste .f1 .f2 .f3 .f4 .f5 | awk '{ 
            printf "%s %s %7.3f\n", $1, $2, $3-($6-$9)+$12+$15 }' > .geo$t
    done
}

annualize() {
    for y in {0..40}; do 
        args=`for m in {0..11}; do printf ".geo%03d " $((12*y+m)); done`
        paste $args | awk '{ a=0; for (i=3;i<=NF;i+=3) a+=$i; print $1" "$2" "a/12 }' > .y$((1979+y))
    done
}

colorize() {
    range=(`sort -nk 3.1 .y* | awk 'NR==1{min=$3} END { print min" "$3 }'`)
    echo ${range[*]}
    for y in {1979..2019}; do awk -vmin=${range[0]} -vmax=${range[1]} 'BEGIN { dlt=max-min }
        {   if ($2 < 191) {$2+=169} else {$2-=191} 
            printf "%s %s %4d\n", $1, $2, 1023*($3-min)/dlt }' .y$y | awk 'BEGIN { n=0
            for (i=255; i>=0; i--) { pal[n] = sprintf("%d 0 255", i); n++ }
            for (i=0; i<=255; i++) { pal[n] = sprintf("0 %d %d", i, 255-i); n++ }
            for (i=0; i<=255; i++) { pal[n] = sprintf("%d 255 0", i); n++ }
            for (i=255; i>=0; i--) { pal[n] = sprintf("255 %d 0", i); n++ }
        } { 
            printf "%s %s %s\n", $2, $1, pal[$3] }
        ' > .c$y
    done
}

scale() {
    rm -f .scale
    range=(`sort -nk 3.1 .y* | awk 'NR==1{min=$3} END { printf "%d %d %d\n", min, $3, $3-min }'`)

    min=${range[0]}; max=${range[1]}; dlt=${range[2]}

    for h in {0..100}; do 
        seq 0 1023 | awk -vh=$h -vmin=$min -vdlt=$dlt 'BEGIN { n=0
                for (i=255; i>=0; i--) { pal[n] = sprintf("%d 0 255", i); n++ }
                for (i=0; i<=255; i++) { pal[n] = sprintf("0 %d %d", i, 255-i); n++ }
                for (i=0; i<=255; i++) { pal[n] = sprintf("%d 255 0", i); n++ }
                for (i=255; i>=0; i--) { pal[n] = sprintf("255 %d 0", i); n++ }
            } { 
            print $1*dlt/1023+min" "h" "pal[$1]
        }' >> .scale 
    done

    echo "set term jpeg size 740,140; set nokey; 
        set title 'Flux (W/m²)
        set xtics 100 out nomirror
        unset ytics; set noborder
        set xrange [$min:$max]; set yrange [0:100]
        rgb(r,g,b) = int(r)*65536 + int(g)*256 + int(b)
        plot '.scale' u 1:2:(rgb(\$3,\$4,\$5)) w dots lc rgb variable lw 1
    " | gnuplot > scale.jpg
}

plot() {
    for y in {1979..2019}; do echo $y >&2; echo "  
        set term jpeg size 740,420; set nokey
        set title '$y'
        set yrange [-180:180]; set xrange [0:720]
        set noborder; unset colorbox 
        unset xtics; unset ytics
        rgb(r,g,b) = int(r)*65536 + int(g)*256 + int(b)
        plot '.c${y}' u (\$1*2):(\$2*2):(rgb(\$3,\$4,\$5)) pt 5 ps 1 lc rgb variable
        " | gnuplot > c$y.jpg
    done
}

animate() {
    convert -loop 0 -delay 50 c*.jpg geoanim.gif
}

clean() { rm -f .geo* .[fyc]* .scale; }

Run it:

$ source geochg.sh
$ require  # Linux Only
$ download
$ extract
$ annualize
$ colorize
$ scale
$ plot
$ animate

Windows users need imagemagick package from Cygwin.