SST

SST_FILES.ersst_file_lists(; path=SST_demo_path)

Build ERSST (Extended Reconstructed SST, monthly, 1854–present) file lists, mirroring file_lists (which targets daily OISST files instead): writes "ersst_whole_file_list.csv" (all months through the latest available) and "ersst_to_get_file_list.csv" (months not yet downloaded) to path.

Source URL pattern: https://www.ncei.noaa.gov/pub/data/cmb/ersst/v5/netcdf/ersst.v5.YYYYMM.nc.

Returns (fil1, fil2), the paths to the whole-list and to-get-list CSVs.

file_lists(path="")

Create file lists and output to csv.

• whole_file_list.csv : all files through today's date
• to_get_file_list.csv : files that remain to download

Sample file names :

url="https://www.ncei.noaa.gov/thredds/dodsC/OisstBase/NetCDF/V2.1/AVHRR/198201/oisst-avhrr-v02r01.19820101.nc"
url="https://www.ncei.noaa.gov/thredds/fileServer/OisstBase/NetCDF/V2.1/AVHRR/198201/oisst-avhrr-v02r01.19820101.nc"

SST_FILES.monthlymean(gdf, m; path0=pwd(), varname="sst")

Compute the mean of varname across all files in group m of grouped file list gdf (e.g. grouped by calendar month), reading each file relative to path0.

Used by SST_processing.monthly_climatology to average all years' daily files for a given calendar month into one climatological monthly field.

SST_FILES.read_map(; variable="anom", file="", file_climatology="")

Read a single day's OISST field from file (falling back to file[1:end-3]*"_preliminary.nc" if file itself doesn't exist — OISST publishes near-real-time files under a _preliminary suffix before the finalized file is available).

variable selects what to return:

• "sst": the raw SST field;
• "anom" (default): the precomputed anomaly field as stored in file;
• "anom_recompute": SST minus the corresponding calendar month's field from file_climatology (read at file's own month, mon_sst, derived from file's date-stamped name) — i.e. an anomaly recomputed against a caller-supplied climatology rather than using the file's own stored anomaly.

Returns the selected 2D field.

test_files(list,ii=[])

Test whether all downloaded files are valid.

list=CSV.read("oisst_whole_file_list.csv",DataFrame)
list_pb=sst_files.test_files(list)
[Downloads.download(r.url,r.fil) for r in eachrow(list[list_pb,:])]

SST_processing.coarse_grain(; datname="oisst", varname="sst", dlon=10.0,
                              path=SST_demo_path, short_demo=false)

Coarse-grain the downloaded OISST daily files (see download_files) onto a dlon-degree grid, writing one CSV per input file plus a single merged CSV of all coarse-grained values.

Reads the file list "$(datname)_whole_file_list.csv" from path (when short_demo is true, only the most recent 10 files); determines the sparse set of non-empty coarse cells once via SST_coarse_grain.indices(list) (reused for every file, since land/ocean geography doesn't change over time). Distributes files evenly across nworkers(): for each file, loads varname from the NetCDF (falling back to the "_preliminary.nc" filename if the primary is missing), computes the coarse-cell area means via SST_coarse_grain.areamean, and writes the result to its own CSV under "$(varname)_lowres_files/$(varname)_lowres_<date>.csv" (any pre-existing output directory for varname is moved aside via mv to a temp path before starting, rather than merged into).

After all files are processed, calls SST_coarse_grain.merge_files to concatenate the per-file CSVs into a single "lowres_oisst_$(varname)_$(dlon).csv".

SST_processing.download_files(; path=SST_demo_path, short_demo=false, verbose=false)

Download the OISST daily NetCDF files listed by SST_FILES.file_lists/ SST_FILES.read_files_list, distributing work across available Julia workers.

If path doesn't exist, it's created. The file list is regenerated via SST_FILES.file_lists(path=path); when short_demo is true, only the most recent 30 files are downloaded (for quick testing). Work is split evenly across nworkers() via @distributed; each missing file is downloaded via Downloads.download, falling back to a "_preliminary.nc"-suffixed URL/filename if the primary download fails (OISST publishes near-real-time files under a _preliminary suffix before the finalized file is available), and silently skipping (with an optional verbose message) if neither is found.

Returns the list of successfully-available local file paths (preferring the finalized file over the preliminary one where both exist), excluding any still-missing entries.

SST_processing.monthly_climatology(; datname="oisst", varname="sst", path=SST_demo_path)

Compute the 1992–2011 monthly climatology (mean SST and mean anomaly, per calendar month) from the downloaded OISST daily files, and write it to a single climatology file via write_SST_climatology.

Reads the file list "$(datname)_whole_file_list.csv" from path, selects the 1992–2011 subset, and groups it by calendar month. For each of "sst" and "anom", computes the across-years mean for each of the 12 calendar months via SST_FILES.monthlymean (distributed across nworkers()), writing each month's mean field via to_monthly_file to a temporary output directory.

Finally combines the 12 per-month files into the single climatology file via write_SST_climatology(output_path, year0, year1, lon, lat) (year0=1992, year1=2011), and returns that file's path.

SST_coarse_grain.areaintegral(arr, i::Int, j::Int, G::NamedTuple, dnl)

Area-weighted sum of arr over the dnl × dnl block of native-resolution cells corresponding to coarse-grid cell (i,j), weighted by the native grid's mask G.msk and cell area G.area. Used by calc_zm to build per-latitude-band area weights.

SST_coarse_grain.areamean(arr, ii, jj, dnl)

Mean of arr over the dnl × dnl block of native-resolution cells corresponding to coarse-grid cell (ii,jj), skipping missing values.

SST_coarse_grain.calc_zm(G::NamedTuple, df, dnl=missing)

Compute a zonal-mean (latitude-band) time series from a coarse-grained, long-format SST table df (as produced by SST_processing.coarse_grain

• SST_coarse_grain.lowres_read, with columns i, j, t, sst for

coarse-grid longitude/latitude indices, time index, and SST value).

G is the coarse grid NamedTuple from grid (lon, lat, msk, area). dnl is the coarse-graining factor in grid cells (e.g. dlon/0.25 for a dlon-degree coarse cell); if not given, defaults to the equivalent of dlon = 10.0.

For each coarse latitude-band index k (from minimum(df.j) to maximum(df.j)), computes the area-weighted mean SST across all longitude cells in that band, for every time step (grouped via groupby(df, :t)). Returns an (nlat, ntime) array arr, with rows for latitude bands outside [minimum(df.j), maximum(df.j)] left as NaN.

cellarea(lon0,lon1,lat0,lat1)

source

As a consequence of a theorem of Archimedes, the area of a cell spanning longitudes l0 to l1 (l1 > l0) and latitudes f0 to f1 (f1 > f0) is

(sin(f1) - sin(f0)) * (l1 - l0) * R^2

where

• l0 and l1 are expressed in radians (not degrees or whatever).
• l1 - l0 is calculated modulo 2*pi (e.g., -179 - 181 = 2 degrees, not -362 degrees).
• R is the authalic Earth radius, almost exactly 6371 km.

grid(fil)

Return (lon=lon,lat=lat,msk=msk,area=area) based on fil.

SST_coarse_grain.indices(list, dlon=10.0)

Determine which coarse-grid cells (at dlon-degree resolution) contain valid (non-NaN) ocean data, using the first file in list as a representative sample.

Returns (i=ii[kk], j=jj[kk], k=kk): the coarse-grid i/j indices of valid cells, and their linear index k into the full coarse grid — used by calc_zm and SST_processing.coarse_grain to avoid computing or storing land/all-NaN cells.

SST_coarse_grain.lowres_index(lon0, lat0, kdf)

Find the index into grouped-keys kdf (coarse-grid (i,j) pairs) whose cell center is nearest to (lon0,lat0).

See also lowres_position (the inverse: index → coordinates).

SST_coarse_grain.lowres_position(ii, jj, kdf)

Convert coarse-grid indices ii,jj (as found in kdf) to their cell- center (longitude, latitude) coordinates, given the coarse resolution implied by kdf's index range (dx = 360/maximum(i)).

See also lowres_index (the inverse: coordinates → index).

lowres_read(;path=SST_demo_path,fil="lowres_oisst_sst_10.0.csv")

Read sst_lowres.csv

lowres_merge(;path=SST_demo_path,variable="sst")

Merge all files found in chosen path.

SST_timeseries.calc(input, list; title="", gdf=nothing)

Compute a full SST time-series diagnostic NamedTuple — raw values, day-of-year climatology, anomaly, and extreme-warm quantile bands — suitable for the Makie extension's by_time/by_year/MHW/ local_and_global plots (via X.options.timeseries).

input is either a raw SST vector, or a DataFrames.GroupKey into gdf (a grouped DataFrame, e.g. grouped by grid cell (i,j)) — in the latter case gdf[input].sst supplies the series. list is the full file/date list (as from SST_FILES.read_files_list), used to align year/month/day with each point in the series and to compute the climatology.

Internally:

1. repeatclim computes the 1992–2011 day-of-year climatology (via clim/gdf_clim) and repeats it across the full series length, giving clim.
2. anom computes sst - climatology, re-centered to the climatology's own median (so anom's scale matches sst, not a zero-centered anomaly).
3. calc_quantile computes, for each day of year, the 10th/90th percentile of the 1992–2011 anomaly (a ±2-day window around each calendar day) — returned as low/high, used to flag extreme warm/cool periods (e.g. the Makie extension's MHW plot).

Returns (sst, clim, anom, title, year, month, day, low, high). title defaults to "SST time series" unless overridden.

SST_plots.MHW(X::SSTdiag)

Plot the SST anomaly time series (X.options.timeseries.sst .- X.options.timeseries.clim), highlighting marine heat wave (MHW) periods in red.

A point is classified as an MHW day when the anomaly exceeds X.options.timeseries.high — the day-of-year 90th percentile threshold computed by SST_timeseries.calc_quantile (see that function for how .high is derived); all other points are blue. period, read via getopt with default (1982,2024), sets the x-axis limits.

SST_plots.TimeLat(X::SSTdiag)

Plot a time-versus-latitude filled-contour anomaly diagram directly from X.options.

Unlike the ECCO extension's TimeLat — which renders a NamedTuple precomputed by a separate ECCO_procs.TimeLat step — this method reads X.options directly, since SST currently has no _procs-equivalent precompute stage. Fields read (via getopt where optional):

• X.options.timeseries — used only for .year's length, to build the time axis
• X.options.zonal_mean — the raw latitude × time array to contour
• period::Tuple = (1982,2024) — x-axis limits
• ylims::Tuple = (-90,90) — y-axis (latitude) limits
• clip_to_range::Bool = true — clip the field to the fixed contour levels (-2.0:0.25:2.0)/5 via to_range! before contouring

Latitude bins are inferred from size(zonal_mean,1), assuming they evenly tile -90:90. The time axis/field are subsampled by a factor of 7 (x[1:7:end], z[1:7:end,:]) for plotting performance.

SST_plots.by_time(X::SSTdiag)

Plot the SST time series X.options.timeseries, optionally overlaid with its seasonal climatology and/or anomaly.

Reads, via getopt (all optional, with the defaults shown):

• X.options.timeseries — a NamedTuple with fields sst, clim, anom, title (as produced by, e.g., SST_timeseries.calc)
• show_anom::Bool = true — overlay timeseries.anom in red
• show_clim::Bool = true — overlay timeseries.clim in orange
• period::Tuple = (1982,2024) — x-axis limits

The time axis assumes timeseries.sst is a daily series starting in 1982 (collect(1:length(ts.sst))/365.25 .+ 1982), independent of period.

SST_plots.by_year(X::SSTdiag)

Plot the SST time series X.options.timeseries.sst as one overlaid line per calendar year, from 1982 through the most recent full year, using a fixed 365-day-per-year assumption (leap days are not separately handled).

Color coding is hardcoded by year:

• 1982–2020: gray
• 2021–2022: blue
• 2023: orange
• 2024 onward: red (linewidth doubled)

These year cutoffs are literal constants in the function body and will need updating in future years to keep highlighting the most recent data; they do not derive from X.options or the current date.

SST_plots.map_base()

Build a base Figure/Axis showing the Blue Marble Next Generation basemap image, for use as a background under subsequent heatmap!/ scatter! calls (see plot_sst_map).

Returns (fig, ax, im).

The image is fetched via MeshArrays.mydatadep("basemap_jpg1"), then reversed, permutedimsd, and circshifted by (1800,0) pixels to align its native orientation and longitude origin with the [0,360) convention used elsewhere in this package (displayed on -0.05 .. 359.95 in x, -89.95 .. 89.95 in y). Axis decorations are hidden.

SST_plots.plot_sst_map(X::SSTdiag)

Plot an SST (or anomaly) map over the map_base basemap, with an optional grid overlay and a highlighted point of interest.

Reads X.options.map_data, a NamedTuple expected to provide: lon, lat, field (the gridded values to show via heatmap!), colormap, colorrange, showgrid::Bool (whether to overlay lowres_scatter's coarse-grain grid-index labels), and lon1/lat1 (coordinates of a single point, marked with both a blue circle and a yellow X — e.g. to indicate the location a companion time-series plot corresponds to).

ERA5_plot.plot_Qnet_cumsum(df, tim, sst)

Plot the normalized cumulative sum of the net surface heat flux anomaly, as a single-panel diagnostic of persistent warming/cooling bias over the record.

Computed as:

z = rnmn(df.qnet, 24)             # 24-hour rolling mean of Qnet
z = cumsum(z .- mean(z))          # cumulative sum of the mean-removed series
z = z ./ sqrt(mean(z.^2))         # normalized to unit RMS ("non-dimensional")

A steadily increasing or decreasing trend indicates a persistent net warming or cooling bias in df.qnet over the plotted period, rather than noise around zero.

ERA5_plot.plot_bulk_formulae(df::DataFrame)
ERA5_plot.plot_bulk_formulae(fil::String)

Plot a stacked column of time series for ["dlw","dsw","hl","hs","qnet"] from bulk-formula output df (one panel per variable, in that order).

The String method is a convenience wrapper: it loads df from file fil via read_bulk_formulae, then delegates to the DataFrame method.

ERA5_plot.plot_surface_balance(df, tim, sst)

Plot a 4-panel surface heat budget summary from bulk-formula output df (e.g. from read_bulk_formulae/surface_balance), a time axis tim (days since Jan. 1), and an SST series sst.

Panels (row, column):

1. (1,1) temperature — df.tmp2m_degC versus sst
2. (1,2) radiative components — 24-hour rolling means (via rnmn) of lw, sw, dlw, dsw, ulw, usw; y-axis fixed to (-300,500) W/m²
3. (2,1) turbulent components & net — rolling means of hl, hs, qnet; y-axis fixed to (-400,400) W/m²
4. (2,2) radiative components & net — rolling means of lw, sw, qnet; y-axis fixed to (-500,300) W/m²

The y-axis ranges in panels 2–4 are fixed constants chosen for visual comparability across components, not derived from df, and may clip series with larger excursions.

E(dtas)

Saturation water vapor pressure via Teten's formula, given temperature dtas (Kelvin). Constants a1=611.21 Pa, a3=17.502, a4=32.19 K, T0=273.16 K.

interpolate_sst(sst, tim)

Linearly interpolate a 365-day daily SST climatology sst (assumed indexed at day-of-year midpoints 0.5:364.5) onto arbitrary time points tim, extrapolating linearly (Line()) beyond the endpoints.

qsat(ps, E)

Saturation specific humidity (equivalent to huss) at surface pressure ps given saturation vapor pressure E (e.g. from E), via (Rdry/Rvap)*E/(ps - (1-Rdry/Rvap)*E), with Rdry=287.0597, Rvap=461.5250 J/(kg·K).

read_bulk_formulae(fil::String)

Demo/sanity-check helper: reads ERA5 variables at the fixed sample point (lon0,lat0) = (205,45) from fil via read_from_nc, computes bulk-formula turbulent fluxes (hl, hs, evap) and net flux qnet against a hardcoded constant SST of 15.0 °C (not a real SST product) — intended for illustrating/testing the bulk-formula pipeline at a single point, not for general-purpose flux computation (see surface_balance for that, which takes sst as an argument).

Returns the augmented DataFrame.

read_from_nc(fil::String, ii, jj)

Read ERA5 reanalysis variables at grid point (ii,jj) from NetCDF file fil, apply unit conversions, and derive specific humidity and wind speed.

Reads 10 ERA5 dataset variables (list_ds, ECMWF short names) and renames/rescales them to list_in (this package's working names), via parallel offset/factor lookup tables indexed by position:

then derives:

• spfh: saturation specific humidity at 2 m, via qsat(pres, E(d2m)) (Teten's formula, see E/qsat);
• wspeed: 10 m wind speed, via wspeed.(u10m, v10m).

Returns a DataFrame with one row per time record in fil, columns list_in plus spfh, wspeed.

read_lonlat(; path_to_data="ERA5_data")

Read the longitude/latitude grid vectors from a fixed sample file, joinpath(path_to_data, "2023/ERA5_2023_01.nc") — i.e. this assumes that specific file exists under path_to_data and that all ERA5 files share the same grid (a reasonable assumption for a fixed reanalysis product, but note the file/path is hardcoded rather than derived from path_to_data generically).

Returns (lon, lat).

using NCDatasets, Climatology
ERA5.read_lonlat()

read_one_year(year,ii,jj)

import Climatology.ERA5: read_lonlat, read_one_year, ij

lon,lat=read_lonlat()

year0=2023
lon0=205; lat0=45;
(ii,jj)=ij(lon0,lat0,lon,lat)

df=read_one_year(year0,ii,jj)

fil="ERA5_lon"*string(lon0)*"_lat"*string(lat0)*"_year"*string(year0)*".csv"
CSV.write(joinpath(tempdir(),fil),df)

using CairoMakie
da=SurfaceFluxDiag((default=true,),df)
plot(da)

read_sample(path=ERA5_OISST_download())

surface_balance(df, sst)

Compute the full air-sea surface heat budget for df given sea surface temperature(s) sst, mutating df in place with the added flux columns and returning it.

Requires df to already have tmp2m_degC, spfh, wspeed (e.g. from read_from_nc) and dlw, dsw (downward long-/shortwave, sign convention: made positive-down via abs.(...) here regardless of input sign).

Computes, via AirSeaFluxes.bulkformulae (turbulent fluxes) plus Stefan-Boltzmann/albedo (radiative fluxes):

• hl, hs, evap: latent heat flux, sensible heat flux, evaporation, from bulkformulae(tmp2m_degC + 273.16, spfh, wspeed, sst) (temperature converted to Kelvin);
• ulw = σ(sst+273.15)^4: upward longwave (Stefan-Boltzmann, σ = 5.670e-8 W·m⁻²·K⁻⁴);
• usw = albedo * |dsw|: upward (reflected) shortwave, albedo = 0.06;
• lw = dlw - ulw, sw = dsw - usw: net longwave/shortwave;
• qnet = hl + hs + lw + sw: total net surface heat flux into the ocean.