Trend in Gross Primary Production

# Zoe Phin, 2021/05/03
# File: prod.sh
# Run: source prod.sh; download; index; plot
# Output: trend info & prod.png

require() { sudo apt-get install -y netpbm gmt; }

sets() {
    for y in {2001..2020}; do
        let n=0
        for w in {0..45}; do
            d=$(date -d "$y-01-01 +$n days" +%Y-%m-%d)
            let n+=8
            echo "wget -O $d.png 'https://gibs.earthdata.nasa.gov/wms/epsg4326/all/wms.cgi?REQUEST=GetMap&SERVICE=WMS&FORMAT=image/png&VERSION=1.1.1&SRS=epsg:4326&BBOX=-180,-90,180,90&TRANSPARENT=TRUE&WIDTH=360&HEIGHT=180&LAYERS=MODIS_Terra_L4_Gross_Primary_Productivity_8Day&TIME=$d'"
        done
    done > sets.sh
}

download() { 
    sets; bash sets.sh; 
    find . -name \*png -type f -size -10k -exec rm {} \;
}

scale() {
    awk 'BEGIN{
        for (i=100;i<=255;i++) printf "%03d 000 000\n", i
        for (i=000;i<=255;i++) printf "255 %03d 000\n", i
        for (i=255;i>=000;i--) printf "%03d 255 000\n", i
        for (i=000;i<=255;i++) printf "000 255 %03d\n", i
        for (i=255;i>=000;i--) printf "000 %03d 255\n", i
    }' | awk '{printf "s/%s/%0.6f/\n", $0, (NR-1)/1180*0.095 }' > replace.sed
    echo "s/000 000 000/999/" >> replace.sed
}

onefile() {
    scale; pngtopnm $1 | pnmtoplainpnm | sed '1,3d;s/  /\n/g' | awk '{
    printf "%03d %03d %03d\n", $1, $2, $3}' | sed -f replace.sed | sed 's/... ... .../999/' | awk '$1!=999{
        l=sprintf("%d",(NR-1)/360)-89.5
        a=6378.137; e=1-6356.752^2/a^2; r=atan2(0,-1)/180; 
        A=(a*r)^2*(1-e)*cos(r*l)/(1-e*sin(r*l)^2)^2
        SA+=A; S+=$1*A
    } END {
        printf "%.6f\n", S/SA
    }'
}

index() {
    for f in $(ls -1 2*.png); do
        echo -n "${f/.png/} "
        onefile $f
    done | tee .csv
    cp .csv bak.csv
}

linear() {
    cat .csv | sed \$d | awk '{ "date +%Y\\ %j -d "$1 | getline t; print t" "$2 }' | awk '
        {printf "%.4f %s\n", $1+$2/365, $3}' | gmt gmtregress | awk '
        NR>1 { printf "%.6f\n", $3 }' | tee .lin | sed -n '1p;$p' | tr '\n' ' ' | awk '{
        printf "%.4f --> %.4f, is %+0.2f%\n", $1, $2, ($2/$1-1)*100 }'
}

plot() { 
    linear; paste -d ' ' .csv .lin > plot.csv
    echo "set term png size 740,620
        set key outside top center horizontal reverse
        set timefmt '%Y-%m-%d'
        set xdata time
        set xtics format '%Y'
        set ytics format '%4.2f'
        set ytics 0.01; set mxtics 2; set mytics 5
        set xrange ['2000-12-01':'2020-02-01']
        set grid xtics mxtics ytics mytics
        plot 'plot.csv' u 1:2 t 'Gross Primary Production (kg C/m²)' w lines lw 2 lc rgb '#00CC00',\
                     '' u 1:3 t 'Linear Regression' w lines lw 3 lc rgb '#005500'		
    " | gnuplot > prod.png 
}