Add nao and sio functions companion of emnso and keeling funs.

src/GMT.jl

@@ -169,7 +169,7 @@ export
 	psconvert, ps2image, ps2raster, psbbox, psview, use_psview, psbasemap, psbasemap!, psclip, psclip!, pscoast, pscoast!, psevents, pshistogram, pshistogram!,
 	psimage, psimage!, pslegend, pslegend!, psmask, psmask!, psrose, psrose!, psscale, psscale!, pssolar, pssolar!,
 	psternary, psternary!, pstext, pstext!, pswiggle, pswiggle!, psxy, psxy!, psxyz, psxyz!, regress, resetGMT, rose,
-	rose!, sample1d, scatter, scatter!, scatter3, scatter3!, solar, solar!, analemma, enso, keeling,
+	rose!, sample1d, scatter, scatter!, scatter3, scatter3!, solar, solar!, analemma, enso, keeling, nao, soi,
 	sunsetrise, spectrum1d, fft1d, fft2d, fft2d!, rgb2hsv, hsv2rgb, sphdistance, sphinterpolate, sphtriangulate, surface, ternary, ternary!, text, text!,
 	text_record, textrepel, annotate, annotate!, trend1d, trend2d, triangulate, gmtsplit,
 	decorated, vector_attrib, wiggle, wiggle!, xyz2grd, gmtbegin, gmtend, gmthelp, subplot, gmtfig, inset, showfig,

src/extras/analemma.jl

@@ -267,3 +267,134 @@ function enso(; data::Bool=false, data0::Bool=false, kwargs...)
 	wiggle(D, track=:faint, ampscale=1.25, figsize=(14,4), R=:tightx, fill=["red+p", "blue+n"], pen=0,
 	       xlabel="Year", title="Oceanic Niño Index", show=true; kwargs...)
 end
+
+# ---------------------------------------------------------------------------------------------------
+"""
+	nao(; data::Bool=false, data0::Bool=false, kwargs...)
+
+Retrieve NAO (North Atlantic Oscillation) index data.
+
+Data is fetched from NOAA (https://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/norm.nao.monthly.b5001.current.ascii),
+the standardized monthly index since 1950. If plotted, the positive phase is shown in red and the negative in blue.
+A plot is generated by default unless `data` or `data0` is set to `true`, case in which the index data
+is returned and no figure is generated.
+
+### Arguments
+- `data::Bool=false`: If `true`, return the NAO data as a [date index] pair in a GMTdataset.
+- `data0::Bool=false`: If `true`, similar to above, but return a 3 column matrix with second column
+   all equal to zero. This is useful for plotting purposes when using the `wiggle` function.
+- `kwargs...`: Additional keyword arguments passed to underlying plotting function.
+
+The default plotting command is:
+```julia
+wiggle(D, track=:faint, ampscale=maximum(abs, D[:,3])/2, figsize=(14,4), R=:tightx, fill=["red+p", "blue+n"], pen=0,
+       xlabel="Year", title="North Atlantic Oscillation Index", show=true; kwargs...)
+```
+where `D` is the dataset returned when `data0=true` (the `ampscale` makes the largest swing 2 cm).
+
+### Returns
+NAO data indices in a GMTdataset or `nothing` depending on the `data` and `data0` flags.
+
+### Examples
+```julia
+nao()			# Plot the NAO index
+```
+"""
+function nao(; data::Bool=false, data0::Bool=false, kwargs...)
+	# The file has one row per month: YEAR MONTH INDEX
+	fname = noaa_climate_download("https://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/norm.nao.monthly.b5001.current.ascii", "nao")
+	M = try
+		gmtread(fname, data=true).data
+	finally
+		rm(fname, force=true)
+	end
+	D = climate_index_ds(M[:,1] .+ (M[:,2] .- 0.5) ./ 12, Float64.(M[:,3]), "NAO", data)
+	(data || data0) && return D
+	wiggle(D, track=:faint, ampscale=climate_ampscale(D), figsize=(14,4), R=:tightx, fill=["red+p", "blue+n"], pen=0,
+	       xlabel="Year", title="North Atlantic Oscillation Index", show=true; kwargs...)
+end
+
+# ---------------------------------------------------------------------------------------------------
+"""
+	soi(; data::Bool=false, data0::Bool=false, kwargs...)
+
+Retrieve SOI (Southern Oscillation Index) data.
+
+Data is NOAA CPC's monthly index, fetched from NOAA PSL (https://psl.noaa.gov/data/correlation/soi.data).
+If plotted, the positive SOI (La Niña) is shown in blue and the negative (El Niño) in red, the same
+colours `enso` gives those phases.
+A plot is generated by default unless `data` or `data0` is set to `true`, case in which the index data
+is returned and no figure is generated.
+
+### Arguments
+- `data::Bool=false`: If `true`, return the SOI data as a [date index] pair in a GMTdataset.
+- `data0::Bool=false`: If `true`, similar to above, but return a 3 column matrix with second column
+   all equal to zero. This is useful for plotting purposes when using the `wiggle` function.
+- `kwargs...`: Additional keyword arguments passed to underlying plotting function.
+
+The default plotting command is:
+```julia
+wiggle(D, track=:faint, ampscale=maximum(abs, D[:,3])/2, figsize=(14,4), R=:tightx, fill=["blue+p", "red+n"], pen=0,
+       xlabel="Year", title="Southern Oscillation Index", show=true; kwargs...)
+```
+where `D` is the dataset returned when `data0=true` (the `ampscale` makes the largest swing 2 cm).
+
+### Returns
+SOI data indices in a GMTdataset or `nothing` depending on the `data` and `data0` flags.
+
+### Examples
+```julia
+soi()			# Plot the SOI
+```
+"""
+function soi(; data::Bool=false, data0::Bool=false, kwargs...)
+	# The file is a "FIRST LAST" years line, then one row per year (YEAR and its 12 monthly values,
+	# -99.99 meaning none), then a few trailer lines. Only rows that are exactly a year are taken.
+	fname = noaa_climate_download("https://psl.noaa.gov/data/correlation/soi.data", "soi")
+	t, v = Float64[], Float64[]
+	try
+		for l in eachline(fname)
+			f = split(l)
+			length(f) != 13 && continue
+			yr = tryparse(Int, f[1])
+			yr === nothing && continue
+			for m = 1:12
+				x = tryparse(Float64, f[m+1])
+				(x === nothing || x <= -99) && continue
+				push!(t, yr + (m - 0.5) / 12);	push!(v, x)
+			end
+		end
+	finally
+		rm(fname, force=true)
+	end
+	isempty(v) && error("The SOI file from NOAA PSL had no monthly values.")
+	D = climate_index_ds(t, v, "SOI", data)
+	(data || data0) && return D
+	wiggle(D, track=:faint, ampscale=climate_ampscale(D), figsize=(14,4), R=:tightx, fill=["blue+p", "red+n"], pen=0,
+	       xlabel="Year", title="Southern Oscillation Index", show=true; kwargs...)
+end
+
+# ---------------------------------------------------------------------------------------------------
+# Download a NOAA climate text file into the user's temp dir. Like in `enso`, curl first (NOAA sites are
+# often very slow with Downloads), Downloads as the fallback. Returns the local file name.
+function noaa_climate_download(url::String, tag::String)::String
+	fname = joinpath(TMPDIR_USR.dir, "GMTjl_" * tag * "_" * TMPDIR_USR.username * TMPDIR_USR.pid_suffix * ".txt")
+	try
+		run(`curl -s --show-error --fail $url -o $fname`)
+	catch
+		Downloads.download(url, fname)
+	end
+	return fname
+end
+
+# A monthly index (T in decimal years, V the index) as `enso` returns its own: [Time index] when DATA,
+# else [Time 0 index] for `wiggle`. The first column is made a Time column.
+function climate_index_ds(T::Vector{Float64}, V::Vector{Float64}, name::String, data::Bool)::GMTdataset
+	D = data ? mat2ds([T V], colnames=["Time", name]) : mat2ds([T zeros(length(T)) V], colnames=["Time", "zero", name])
+	setdecyear_time!(D)
+	return D
+end
+
+# The wiggle's amplitude scale (index units per cm) that makes the largest swing of D's index 2 cm, half
+# the 4 cm tall figure, so no peak is cut by the frame. `enso`'s fixed 1.25 is this for the ONI (~2.6 peak).
+climate_ampscale(D::GMTdataset)::Float64 = round(maximum(abs, view(D.data, :, size(D, 2))) / 2, digits=2)

test/test_misc.jl

@@ -436,5 +436,13 @@
 	analemma(show=false)
 	keeling(show=false)
 	enso(show=false)
+	D = nao(data=true)
+	@test D.colnames == ["Time", "NAO"] && size(D, 1) > 900 && D[1,1] ≈ 1950 + 0.5/12
+	@test size(nao(data0=true), 2) == 3
+	nao(show=false)
+	D = soi(data=true)
+	@test D.colnames == ["Time", "SOI"] && size(D, 1) > 800 && all(D[:,2] .> -99)
+	@test size(soi(data0=true), 2) == 3
+	soi(show=false)
 
 end