(45) Trend modeling

In this example show how the module trend1d is used to fit the CO2 data set collected from the top of Mauna Loa. This yields the famous Keeling curve.

#!/usr/bin/env bash
#               GMT EXAMPLE 45
#
# Purpose:      Illustrate use of trend1d mixed models
# GMT modules:  text, plot, trend1d
# Unix progs:   rm
#

gmt begin ex45
	gmt set FONT_HEADING 24p
	gmt subplot begin 4x1 -Fs15c/4.8c -R1958/2016/310/410 -Bxaf -Byaf+u" ppm" -BWSne+gazure1 -T"The Keeling Curve [CO@-2@- at Mauna Loa]"

		# Basic LS line y = a + bx
		gmt subplot set 3
		gmt trend1d -Fxm @MaunaLoa_CO2.txt -Np1 > model.txt
		gmt plot -Sc0.05c -Gred @MaunaLoa_CO2.txt
		gmt plot -W0.5p,blue model.txt
		gmt text -F+f12p+cTL+t"m@-2@-(t) = a + b\267t" -Dj8p -Glightyellow

		# Basic LS line y = a + bx + cx^2
		gmt subplot set 2
		gmt trend1d -Fxm @MaunaLoa_CO2.txt -Np2 > model.txt
		gmt plot -Sc0.05c -Gred @MaunaLoa_CO2.txt
		gmt plot -W0.5p,blue model.txt
		gmt text -F+f12p+cTL+t"m@-3@-(t) = a + b\267t + c\267t@+2@+" -Dj8p -Glightyellow

		# Basic LS line y = a + bx + cx^2 + seasonal change
		gmt subplot set 1
		gmt trend1d -Fxmr @MaunaLoa_CO2.txt -Np2,f1+o1958+l1 > model.txt
		gmt plot -Sc0.05c -Gred @MaunaLoa_CO2.txt
		gmt plot -W0.25p,blue model.txt
		gmt text -F+f12p+cTL+t"m@-5@-(t) = a + b\267t + c\267t@+2@+ + d\267cos(2@~p@~t) + e\267sin(2@~p@~t)" -Dj8p -Glightyellow

		# Plot residuals of last model
		gmt subplot set 0
		gmt plot -R1958/2016/-4/4 -Byafg10+u" ppm" -Sc0.05c -Gred model.txt -i0,2
		gmt text -F+f12p+cTL+t"@~e@~(t) = y(t) - m@-5@-(t)" -Dj8p -Glightyellow

	gmt subplot end
	rm -f model.txt
gmt end show