diff --git a/src/GMT.jl b/src/GMT.jl index 34dac3b13..87760b098 100644 --- a/src/GMT.jl +++ b/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, diff --git a/src/extras/analemma.jl b/src/extras/analemma.jl index 84b78b2bd..2a4f45a5c 100644 --- a/src/extras/analemma.jl +++ b/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) diff --git a/test/test_misc.jl b/test/test_misc.jl index 2c710d28d..f87239b67 100644 --- a/test/test_misc.jl +++ b/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