|
| 1 | +""" |
| 2 | + analemma(; lon=0, lat=0, hour=12, year=2024, lonlat=false, data=false, cmap=:turbo, kwargs...) |
| 3 | +
|
| 4 | +Plot the analemma. |
| 5 | +
|
| 6 | +The analemma is the figure-8 pattern traced by the Sun's position in the sky |
| 7 | +when observed from the same location at the same mean solar time throughout the year. |
| 8 | +
|
| 9 | +### Arguments |
| 10 | +- `lon` : Observer longitude in degrees. |
| 11 | +- `lat` : Observer latitude in degrees. If both `lon` and `lat` are equal to zero (the default) |
| 12 | + we compute your approximate location using web service based your IP address |
| 13 | + (this adds a delay to the calculations). |
| 14 | +- `year` : Year for calculation (default: 2026). This is of minor importance but things slowly change with time. |
| 15 | +- `hour` : Local mean solar time hour (0-24, default: 12 = noon) |
| 16 | +- `cmap` : Colormap for day-of-year coloring (default: :turbo) |
| 17 | +- `lonlat=false`: If `true`, plot longitude vs latitude; if `false` (default), plot azimuth vs elevation. |
| 18 | +- `data=false`: If `true`, return the CO2 data in a GMTdataset. |
| 19 | +- Additional kwargs are passed to `plot` |
| 20 | +
|
| 21 | +### Example |
| 22 | +```julia |
| 23 | +analemma(lon=-9, lat=38.7) # Noon analemma in Lisbon |
| 24 | +analemma(lat=-23.5, hour=9) # 9 AM analemma in tropics |
| 25 | +
|
| 26 | +D = analemma(lon=-9, lat=38.7, lonlat=true, data=true) # Get analemma in lon,lat as a GMTdataset |
| 27 | +``` |
| 28 | +""" |
| 29 | + |
| 30 | +function analemma(; lon::Real=0, lat::Real=0, hour::Real=12, year::Int=2026, |
| 31 | + lonlat=false, data=false, cmap=:turbo, kwargs...) |
| 32 | + d = KW(kwargs) |
| 33 | + if (lon == 0 && lat == 0) lon, lat = get_my_lonlat() end |
| 34 | + _analemma(Float64(lon), Float64(lat), year, Float64(hour), cmap, lonlat==1, data==1, d) |
| 35 | +end |
| 36 | + |
| 37 | +function _analemma(lon, lat, year, hour, cmap, lonlat::Bool, data::Bool, d) |
| 38 | + |
| 39 | + TZ = round(Int, 90 / 15) # Approximate time zone from longitude |
| 40 | + n_days = Dates.isleapyear(year) ? 366 : 365 |
| 41 | + ana = Matrix{Float64}(undef, n_days, 4) |
| 42 | + k = lonlat ? (1,2,3,4) : (3,4,1,2) # lonlat: lon, lat, az, el ; else: az, el, lon, lat |
| 43 | + cnames = lonlat ? ["Lon", "Lat", "Azimuth", "Elevation"] : ["Azimuth", "Elevation", "Lon", "Lat"] |
| 44 | + |
| 45 | + hhmm = @sprintf("T%02d:%02d:00", floor(Int, hour), round(Int, getdecimal(hour) * 60)) |
| 46 | + |
| 47 | + for day in 1:n_days |
| 48 | + date = Date(year, 1, 1) + Dates.Day(day - 1) |
| 49 | + datetime_str = Dates.format(date, "yyyy-mm-dd") * hhmm |
| 50 | + |
| 51 | + result = gmt(@sprintf("solar -I%g/%g+d%s+z%d -C", lon, lat, datetime_str, TZ)) |
| 52 | + ana[day, 1], ana[day, 2], ana[day, 3], ana[day, 4] = result.data[k[1]], result.data[k[2]], result.data[k[3]], result.data[k[4]] |
| 53 | + end |
| 54 | + |
| 55 | + D = mat2ds(ana) |
| 56 | + data && (D.colnames = cnames; return D) |
| 57 | + |
| 58 | + do_show = ((val = find_in_dict(d, [:show])[1]) === nothing) ? true : false # Default is to show |
| 59 | + C = makecpt(cmap=cmap, range=(1, n_days)) |
| 60 | + plot(D, marker="c", markersize="4p", cmap=C, zcolor=collect(1:n_days), colorbar=(xlabel="Day of year",), |
| 61 | + xlabel= lonlat ? "Longitude" : "Azimuth", ylabel= lonlat ? "Latitude" : "Elevation", |
| 62 | + title=@sprintf("Analemma %02d:00 lat=%.1f", floor(Int, hour), lat), show=do_show, d...) |
| 63 | +end |
| 64 | + |
| 65 | +# --------------------------------------------------------------------------------------------------- |
| 66 | +""" |
| 67 | + sunsetrise(; lon=0, lat=0, year=2026, TZ::Int=50, raise=false, both=false, data=false; kwargs...) |
| 68 | +
|
| 69 | +Plot sunrise and sunset times throughout the year for a given location. |
| 70 | +
|
| 71 | +Uses GMT's `solar` module for calculations. |
| 72 | +
|
| 73 | +### Arguments |
| 74 | +- `lon` : Observer longitude in degrees. |
| 75 | +- `lat` : Observer latitude in degrees. If both `lon` and `lat` are equal to zero (the default) |
| 76 | + we compute your approximate location using web service based your IP address |
| 77 | + (this adds a delay to the calculations). |
| 78 | +- `year` : Year for calculation (default: 2026). This is of minor importance but things slowly change with time. |
| 79 | +- `TZ` : Time zone offset in hours. By default (when the default value of 50 stands) we compute |
| 80 | + it from longitude but it doesn't take into account daylight saving time. |
| 81 | +- `raise=false`: If `true`, plot sunrise times; if `false`, plot sunset times. |
| 82 | +- `both=false`: If `true`, plot both sunrise and sunset times. |
| 83 | +- `data=false`: If `true`, return the sunset or sunrise data (depending on `rise`) |
| 84 | + or both if `both=true` in a GMTdataset. |
| 85 | +- Additional kwargs are passed to `plot` |
| 86 | +
|
| 87 | +### Returns |
| 88 | +If `data=true` returns a GMTdataset if `both` is not set (false) or a tuple of GMTdatasets with |
| 89 | +sunrise and sunset data if `both=true`. Returns `nothing` if a plot is made. |
| 90 | +
|
| 91 | +Example |
| 92 | +------- |
| 93 | +```julia |
| 94 | + sunsetrise(lat=38.7, lon=-9) # Lisbon |
| 95 | + sunsetrise(lat=60) # High latitude with long summer days |
| 96 | +
|
| 97 | + Dsrise, Dsset = sunsetrise(lat=38.7, lon=-9, both=true, data=true) # Get sunrise/set data |
| 98 | +``` |
| 99 | +""" |
| 100 | +function sunsetrise(; lon=0.0, lat=0.0, year::Int=2026, TZ::Int=50, raise=false, both=false, |
| 101 | + data::Bool=false, kwargs...) |
| 102 | + d = KW(kwargs) |
| 103 | + _TZ = (TZ == 50) ? round(Int, (datetime2unix(now()) - datetime2unix(now(UTC))) / 3600) : TZ |
| 104 | + if (lon == 0 && lat == 0) lon, lat = get_my_lonlat() end |
| 105 | + _sunsetrise(Float64(lon), Float64(lat), year, _TZ, raise==1, both==1, data==1, d) |
| 106 | +end |
| 107 | +function _sunsetrise(lon, lat, year::Int, TZ::Int, raise::Bool, both::Bool, data::Bool, d) |
| 108 | + |
| 109 | + n_days = Dates.isleapyear(year) ? 366 : 365 |
| 110 | + sunrise = Matrix{Float64}(undef, n_days, 2) |
| 111 | + both && (sunset = Matrix{Float64}(undef, n_days, 2)) |
| 112 | + sun = sunrise # For raise or set |
| 113 | + ind = raise ? 5 : 6 # 5 for raise, 6 for set |
| 114 | + |
| 115 | + for day in 1:n_days |
| 116 | + date = Date(year, 1, 1) + Dates.Day(day - 1) |
| 117 | + datetime_str = Dates.format(date, "yyyy-mm-dd") |
| 118 | + |
| 119 | + # Use solar: columns are lon, lat, az, el, sunrise, sunset, noon, duration, ... |
| 120 | + # Values are in fraction of day, multiply by 24 to get hours |
| 121 | + result = gmt(@sprintf("solar -I%g/%g+d%s+z%d -C", lon, lat, datetime_str, TZ)) |
| 122 | + ydec = datetime2unix(yeardecimal(year + (day - 0.5) / n_days)) |
| 123 | + both ? (sunrise[day, 1] = ydec; sunrise[day, 2] = result[1,5] * 24; |
| 124 | + sunset[day, 1] = ydec; sunset[day, 2] = result[1,6] * 24) : |
| 125 | + (sun[day, 1] = ydec; sun[day, 2] = result[1,ind] * 24) |
| 126 | + end |
| 127 | + |
| 128 | + doy = dayofyear(today()) |
| 129 | + |
| 130 | + both ? (Dsr = mat2ds(sunrise); settimecol!(Dsr, col=1); Dss = mat2ds(sunset); settimecol!(Dss, col=1)) : |
| 131 | + (Dsun = mat2ds(sun); settimecol!(Dsun, col=1)) |
| 132 | + |
| 133 | + data && return both ? (Dsr, Dss) : Dsun |
| 134 | + |
| 135 | + title = @sprintf("%s lon=%.2f lat=%.2f", raise ? "Sunrise" : both ? "Sunrise/Sunset" : "Sunset", lon, lat) |
| 136 | + y_label = "Hour (UTC $(TZ))" |
| 137 | + xaxis_nt = (axes=:Sen, annot=1, annot_unit=:month, ticks=7, ticks_unit=:day_date) |
| 138 | + yaxis_nt = (annot=15, annot_unit=:minute2, ticks=5, ticks_unit=:minute2, label=y_label) |
| 139 | + par = (FORMAT_DATE_MAP="o", FORMAT_TIME_PRIMARY_MAP="abbreviated") |
| 140 | + |
| 141 | + do_show = ((val = find_in_dict(d, [:show])[1]) === nothing) ? true : false # Default is to show |
| 142 | + fmt::String = ((val = find_in_dict(d, [:fmt])[1]) !== nothing) ? arg2str(val)::String : FMT[]::String |
| 143 | + savefig = ((val = find_in_dict(d, [:savefig :figname :name])[1]) !== nothing) ? arg2str(val)::String : nothing |
| 144 | + opt_R = ((val = find_in_dict(d, [:R :region :limits])[1]) !== nothing) ? val : "tightx" |
| 145 | + both ? plotyy(Dsr, Dss, yaxis=yaxis_nt, title=title, conf=par, R=opt_R, lw=1; d...) : |
| 146 | + plot(Dsun, xaxis=xaxis_nt, yaxis=yaxis_nt, title=title, lc="#0072BD", lw="1p", conf=par, R=opt_R; d...) |
| 147 | + |
| 148 | + use_back = (CTRL.limits[7] == 0.0 && CTRL.limits[8] == 0.0) # Only used if -Rtight |
| 149 | + back_lims = CTRL.limits[1:4] |
| 150 | + both ? plot!([Dsr[doy:doy,1:2]; Dss[doy:doy,1:2]], marker=:circle, mc=:yellow, ms="6p", mec=:black, fmt=fmt, name=savefig, show=do_show) : |
| 151 | + plot!(Dsun[doy:doy,1:2], marker=:circle, mc=:yellow, ms="6p", mec=:black) |
| 152 | + |
| 153 | + if (!both) |
| 154 | + lims = use_back ? back_lims : CTRL.limits[7:10] |
| 155 | + opt_R=@sprintf("%f/%f/%ft/%ft", lims[1], lims[2], lims[3]/24, lims[4]/24) |
| 156 | + basemap!(frame=(axes=:W, annot="15m", ticks="5m", label=y_label), axis2=(annot=1, annot_unit=:hour), R=opt_R, name=savefig, |
| 157 | + fmt=fmt, conf=(FORMAT_CLOCK_MAP="-hham", FONT_ANNOT_PRIMARY="+9p", TIME_UNIT="d"), show=do_show) |
| 158 | + end |
| 159 | +end |
| 160 | + |
| 161 | +# --------------------------------------------------------------------------------------------------- |
| 162 | +""" |
| 163 | + keeling(; data::Bool=false, kwargs...) |
| 164 | +
|
| 165 | +Plot the Keeling Curve - atmospheric CO2 concentration measured at Mauna Loa since 1958. |
| 166 | +
|
| 167 | +Data is fetched from NOAA (https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.txt). |
| 168 | +
|
| 169 | +### Arguments |
| 170 | +- `data::Bool=false`: If `true`, return the CO2 data in a GMTdataset. |
| 171 | +- Additional kwargs are passed to `plot` |
| 172 | +
|
| 173 | +### Returns |
| 174 | +A GMTdataset of CO2 data if `data=true`, or `nothing` if a plot is made. |
| 175 | +
|
| 176 | +### Examples |
| 177 | +```julia |
| 178 | +D = keeling(data=true) # Get CO2 data as GMTdataset |
| 179 | +
|
| 180 | +keeling(lw=1, lc=:darkgreen) # Plot with custom line width and color |
| 181 | +``` |
| 182 | +""" |
| 183 | +function keeling(; data::Bool=false, kwargs...) |
| 184 | + |
| 185 | + opt_i = data ? "2,3,4,5,6,7" : "2,3" |
| 186 | + D = gmtread("https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.txt", i=opt_i) |
| 187 | + setdecyear_time!(D) # First column is decimal year, make a Time column |
| 188 | + if (data) |
| 189 | + D.colnames[2:end] = ["monthly_average", "de-seasonalized", "#days", "st.dev_of_days", "unc.of_mon_mean"] |
| 190 | + return D |
| 191 | + end |
| 192 | + |
| 193 | + plot(D, lw=0.75, lc=:red, xlabel="Year", ylabel="CO@-2@- (ppm)", title="Keeling Curve - Mauna Loa CO@-2@-", |
| 194 | + show=true; kwargs...) |
| 195 | +end |
| 196 | + |
| 197 | +# --------------------------------------------------------------------------------------------------- |
| 198 | +""" |
| 199 | + enso(; data::Bool=false, data0::Bool=false, kwargs...) |
| 200 | +
|
| 201 | +Retrieve ENSO (El Niño-Southern Oscillation) data. |
| 202 | +
|
| 203 | +Data is fetched from NOAA (https://www.cpc.ncep.noaa.gov/data/indices/oni.ascii.txt). If plotted, |
| 204 | +El Niño events (positive) shown in red, La Niña (negative) in blue. |
| 205 | +A plot is generated by default unless `data` or `data0` is set to `true`, case in which the index data |
| 206 | +is returned and no figure is generated. |
| 207 | +
|
| 208 | +### Arguments |
| 209 | +- `data::Bool=false`: If `true`, return the computed ENSO data as a [date index] pair in a GMTdataset. |
| 210 | +- `data0::Bool=false`: If `true`, similar to above, but return a 3 column matrix with second column |
| 211 | + all equal to zero. This is useful for plotting purposes when using the `wiggle` function. |
| 212 | +- `kwargs...`: Additional keyword arguments passed to underlying plotting function. |
| 213 | +
|
| 214 | +Note, the plot is created with a figure size of (14,4), with x-axis labeled "Year" and title |
| 215 | +"Oceanic Niño Index", but this can be overwritten via the `xlabel` and `title` options. The option |
| 216 | +`data0` returns a dataset with a zero middle column useful for plotting with `wiggle`. The default |
| 217 | +plotting command is: |
| 218 | +```julia |
| 219 | +wiggle(D, track=:faint, ampscale=1.25, figsize=(14,4), R=:tightx, fill=["red+p", "blue+n"], pen=0, |
| 220 | + xlabel="Year", title="Oceanic Niño Index", show=true; kwargs...) |
| 221 | +``` |
| 222 | +where `D` is the dataset returned when `data0=true`. You can use this to customize a new plot further. |
| 223 | +
|
| 224 | +And, it seems that the NOOA site sometimes is quite slow to respond, so be patient! |
| 225 | +
|
| 226 | +### Returns |
| 227 | +ENSO data indices in a GMTdataset or `nothing` depending on the `data` and `data0` flags. |
| 228 | +
|
| 229 | +### Examples |
| 230 | +```julia |
| 231 | +enso() # Plot ENSO index |
| 232 | +``` |
| 233 | +""" |
| 234 | +function enso(; data::Bool=false, data0::Bool=false, kwargs...) |
| 235 | + |
| 236 | + # Fetch data |
| 237 | + resp = Downloads.download("https://www.cpc.ncep.noaa.gov/data/indices/oni.ascii.txt") |
| 238 | + lines = readlines(resp) |
| 239 | + |
| 240 | + year, vals = Float64[], Float64[] |
| 241 | + |
| 242 | + # ONI format: SEAS YR TOTAL ANOM |
| 243 | + # e.g.: DJF 1950 24.72 -1.53 |
| 244 | + month_map = Dict("DJF"=>1, "JFM"=>2, "FMA"=>3, "MAM"=>4, "AMJ"=>5, "MJJ"=>6, |
| 245 | + "JJA"=>7, "JAS"=>8, "ASO"=>9, "SON"=>10, "OND"=>11, "NDJ"=>12) |
| 246 | + |
| 247 | + for line in lines |
| 248 | + contains(line, "SEAS") && continue # Skip header |
| 249 | + isempty(strip(line)) && continue |
| 250 | + parts = split(line) |
| 251 | + length(parts) < 4 && continue |
| 252 | + |
| 253 | + try |
| 254 | + season = String(parts[1]) |
| 255 | + yr = parse(Float64, parts[2]) |
| 256 | + anom = parse(Float64, parts[4]) |
| 257 | + |
| 258 | + # Convert to decimal year |
| 259 | + mon = get(month_map, season, 1) |
| 260 | + dec_year = yr + (mon - 0.5) / 12 |
| 261 | + |
| 262 | + push!(year, dec_year) |
| 263 | + push!(vals, anom) |
| 264 | + catch |
| 265 | + continue |
| 266 | + end |
| 267 | + end |
| 268 | + rm(resp) |
| 269 | + |
| 270 | + D = data ? mat2ds([year vals]) : mat2ds([year zeros(length(year)) vals]) |
| 271 | + setdecyear_time!(D) # First column is decimal year, make a Time column |
| 272 | + D.colnames[data ? 2 : 3] = "ONI" # Will be wrong for plotting but in that case we don't care |
| 273 | + |
| 274 | + (data || data0) && return D |
| 275 | + |
| 276 | + wiggle(D, track=:faint, ampscale=1.25, figsize=(14,4), R=:tightx, fill=["red+p", "blue+n"], pen=0, |
| 277 | + xlabel="Year", title="Oceanic Niño Index", show=true; kwargs...) |
| 278 | +end |
0 commit comments