Regarding the “ice-coverage” plot I mentioned adding …
Do you have a more clever use of grdmath ? Creating a 2-D grid for 1 value doesn’t seem efficient.
projection_north='-Js0/90/77/1:93300000'
region_north='-R0/360/50/90'
for ff in $(ls -v data/*nc);
do
dd=$(ls ${ff} | awk -F[_.] '{print $2}')
gmt grdproject $projection_north $region_north ${ff}?ice_conc -Gdata_bis.nc -C -I -Fk
gmt grdmath data_bis.nc 1 GE SUM = tmp.nc
gmt math -Q tmp.nc?z[0] = | awk -v day=${dd} 'NR==1 {print day,$0}' >> x.txt
done
gmt begin ice_cov png
gmt plot -R2016-01-01T/2016-12-31T/0/150000 -JX10cT/10c x.txt -Wred --FORMAT_DATE_IN=yyyymmdd -B+t"Sea-ice coverage" -l"North Pole" -Bxafg -Byafg+l"cell count"
gmt end show
(Also, converting cell count to km^2 could be a nice thing to do)