Plot -w*flag* bug?

Q&A
1

Hi there,

I have a time-series I’d like to plot.
I start by generating the points (with gmt math, see below) between 1 (January 1st) and 365 (December 31st).

  1. Can someone check where I have this “extra operand” left on the stack ?
  2. For some reason if I specify x to be “time” via -JX10ct/10c, I got date format HH:MM instead of days despite :
    • specifying --TIME_UNIT=d
    • specifying --DATE_FORMAT_OUT=o (for abbreviated month)
  3. Is it possible to shift the x-frame to 0 (middle of the plot) ?

Thanks !

gmt math -T1/365/1 2 PI MUL T 1 SUB MUL 365 DIV STO@r \
0.006818 \
@r COS 0.39912 MUL SUB \
@r SIN 0.0757 MUL ADD \
2 @r MUL COS 0.006758 MUL SUB \
2 @r MUL SIN 0.000907 MUL ADD \
3 @r MUL COS 0.002697 MUL SUB \
3 @r MUL SIN 0.00148 MUL ADD \
180 PI DIV MUL = data.txt

gmt begin test png
    gmt plot data.txt -R-4/9/-25/25 -JX10ct/10c --TIME_UNIT=d -wa
    gmt basemap -Bxaf -Byaf+u@. -BWS+t"Curve" --MAP_FRAME_TYPE=graph --MAP_VECTOR_SHAPE=0.5 --FORMAT_DATE_OUT=o
gmt end show
2

Obviously I could simply create a vector -T-101/264/1 and avoid -w … but where’s the fun in that ?

gmt math -T-99/266/1 2 PI MUL T 1 SUB MUL 365 DIV STO@r \
0.006818 \
@r COS 0.39912 MUL SUB \
@r SIN 0.0757 MUL ADD \
2 @r MUL COS 0.006758 MUL SUB \
2 @r MUL SIN 0.000907 MUL ADD \
3 @r MUL COS 0.002697 MUL SUB \
3 @r MUL SIN 0.00148 MUL ADD \
180 PI DIV MUL = data.txt

gmt begin test png
    gmt plot data.txt -R-99/266/-25/25 -JX10ct/10c -BW+t"Curve"  -Byaf+u@. --MAP_FRAME_TYPE=graph --TIME_UNIT=d 
    gmt basemap -BS -Bxa1Of1O --FORMAT_DATE_MAP=o --FORMAT_TIME_PRIMARY_MAP=a --MAP_FRAME_TYPE=graph --TIME_UNIT=d -Y5c
gmt end show

test
test1374×1292 22.4 KB

Still … if that could be achieved easily with -wflag I’m interested …

3

I think STO@r does not pop anything off the stack, it just saves that shorthand, so you still have what you call @r sitting on the stack. If you add POP after STO@r it fixes it.

I think -wa (as well as the other time wrappers) expect absolute time but you are giving relative time.

Very nice though!

4

Also I checked the gmt math documentation and it does mention this:

Note that STO and CLR leave the stack unchanged

5

I haven’t tried, but do you think using -ft in plot would help ? Or do I simply use absolute time in math?

6

@pwessel
One last bit … sorry it’s for my sanity … Is there an easy way to get rid of the first folded “September” ?
I tried to play with -Ba1O+phase but it was not conclusive…

#!/bin/bash
# basemap, begin, coast, grdimage, grdmath, info, inset, legend, makecpt, math, plot, text

gmt math -foT --TIME_UNIT=d -T0/365/1 2 PI MUL T 1 SUB MUL 365 DIV STO@r POP \
                                0.006818 \
                                @r COS 0.39912 MUL SUB \
                                @r SIN 0.0757 MUL ADD \
                                2 @r MUL COS 0.006758 MUL SUB \
                                2 @r MUL SIN 0.000907 MUL ADD \
                                3 @r MUL COS 0.002697 MUL SUB \
                                3 @r MUL SIN 0.00148 MUL ADD \
                                180 PI DIV MUL = declination.txt

active_day="2012-12-21T"
active_day_rel=$(gmt math -Q $active_day $(gmt math -Q -fT $active_day --FORMAT_DATE_OUT=yyyy --TIME_UNIT=y --FORMAT_CLOCK_OUT=- =)-01-01T SUB --TIME_UNIT=d =)
active_line=$(gmt math -Q $active_day_rel 1 ADD =)

wrap_ref="1970-09-23T"
wrap_shift=$(gmt math -Q $wrap_ref 1970-12-31T SUB --TIME_UNIT=d 365 DIV =)
wrap_close=$(gmt math -Q $wrap_shift 1 ADD =)


wrap_day=$(gmt math -Q $active_day_rel 365 DIV STO@r $wrap_close GT 1 MUL STO@r2 POP @r @r2 SUB =)

delta=$(echo $(awk  -v var=$active_line 'NR==var {print $2}' declination.txt))
delta_max=$(gmt info declination.txt -C -o3)

gmt begin test png

    gmt plot declination.txt -R$wrap_shift/$wrap_close/-25/25 -wy -JX10c -BW  -Byaf+u@. --MAP_FRAME_TYPE=graph
    echo $wrap_day $delta | gmt plot -Sc0.1c -Gred -Wfaint,black
    gmt basemap -R$wrap_ref/$(gmt math -Q -fT $wrap_ref 365 ADD --TIME_UNIT=d =)/-25/25 -BS -Bxa1Of --FORMAT_DATE_MAP=o --FORMAT_TIME_PRIMARY_MAP=a --MAP_FRAME_TYPE=graph -Y5c 
    gmt text -Dj0.2c -N -F+f+j -Y-5c <<- EOF
		1971-06-21T 0 12p,Times-Roman BC Boreal summer
		1970-12-31T 0 12p,Times-Roman BC Austral summer
		$wrap_ref 30 12p,Symbol RT d
	EOF

    gmt inset begin -DjBR+w3.5c+o-3c/0.25c
        
        gmt grdmath -Rg -I0.5 X SIND = illu.nc
        gmt makecpt -Cgray -T-1/0 -I > illu.cpt
		
        gmt coast -Rg -JG0/0/3c+w$(gmt math -Q $delta_max -1 MUL =) -Dc -Ggray -A5000 -Bg -Wfaint -X0.25 -Y0.25
        gmt grdimage -JG0/0/3c+w$(gmt math -Q $delta $delta_max ADD -1 MUL =) illu.nc -C -t25
        gmt plot -JG0/0/3.5c+w$(gmt math -Q $delta_max -1 MUL =) -Wthin,black -Am -X-0.25 -Y-0.25 <<- EOF
			0 -89
			0 89
		EOF
        gmt plot -Sv0.15c+ea+h0+a60+gred -Wthicker,red -N <<- EOF
			-90 $delta $(gmt math -Q $delta $delta_max ADD -1 MUL =) 0.3c 
		EOF

    gmt inset end

    gmt legend -DjTR+o-3c/-1c -F <<- EOF
		H 14p,Helvetica-Bold Solar declination
		D 5p 1p
		S 0.5c c 0.1c red faint 1c Observed day
		S 0.5c v0.15c+h0+ea+a60 0.3c red thicker,red 1c Solar equatorial plan
	EOF
gmt end show
7

I tried to add -Bx+el but it has no effect whatsoever …

(FYI : the option +el works on minimal test)
…

Ok it’s just not compatible with dates apparently