Select, group and reduce data sets with CommonDataModel
As an example we use the NOAA OISST v2 dataset available at: https://psl.noaa.gov/thredds/fileServer/Datasets/noaa.oisst.v2.highres/sst.day.mean.2023.nc
using CommonDataModel: @select, @CF_str, @groupby, groupby
using Dates
using Downloads: download
using CairoMakie # use GLMakie for interactive plots
#using GLMakie
using IntervalSets
using NCDatasets
using Statistics
Some helper functions for plotting with Makie for plotting maps and timeseries.
function nicemaps(v; timeindex = 1, lon = v["lon"][:], lat = v["lat"][:], title = nothing)
fig, ax, hm = heatmap(lon,lat,v[:,:,timeindex],aspect_ratio = 1/cosd(mean(lat)))
if !isnothing(title)
ax.title[] = title
end
Colorbar(fig[1,2],hm)
return fig
end
function nicetimeseries(v::AbstractVector; time = v["time"][:], title = nothing,
timefmt = "yyyy-mm-dd")
fig, ax, ln = lines(Dates.value.(time),v[:])
if !isnothing(title)
ax.title[] = title
end
timeticks = DateTime.(2023,1:2:12,1)
ax.xticks = Dates.value.(timeticks),Dates.format.(timeticks,timefmt)
ax.xticklabelrotation = π/4
return fig
end
nicetimeseries (generic function with 1 method)
Download the data file (unless already downloaded) The dataset contains the variables lon
(longitude), lat
(latitude), time
and sst
(sea surface temperature).
url = "https://psl.noaa.gov/thredds/fileServer/Datasets/noaa.oisst.v2.highres/sst.day.mean.2023.nc"
fname = "sst.day.mean.2023.nc"
if !isfile(fname)
download(url,fname)
end
"sst.day.mean.2023.nc"
Select a subset by values from coordinate variables
Open the dataset and access the variable sst
. Select the domain from 300°E to 360°E and 0°N to 46°N and the closest time instance to the 1 April 2023.
ds = NCDataset(fname)
ncsst = ds["sst"]
ncsst2 = @select(ncsst,300 <= lon <= 360 && 0 <= lat <= 46 && time ≈ DateTime(2023,4,1))
nicemaps(ncsst2, title = "NA SST")
One can also select a subset by indices and by values. Here we load the first time instance and the domain from 300°E to 360°E and 0°N to 46°N
ncsst_first = view(ncsst,:,:,1)
ncsst_na = @select(ncsst_first,300 <= lon <= 360 && 0 <= lat <= 46)
nicemaps(ncsst_na)
Select all data from a given month
ncsst_march = @select(ncsst,Dates.month(time) == 3)
nicemaps(ncsst_march, timeindex = 1, title = "SST 1st March 2023")
Select multiple months using e.g. an interval (from IntervalSets). ∈ can be typed by writing \in
directly followed by the TAB key.
ncsst_march_april = @select(ncsst,Dates.month(time) ∈ 3..4)
nicemaps(ncsst_march_april)
Use julia functions to extract data; here the abs
function polar regions north of 60°N and south of 60°S
ncsst_polar = @select(ncsst,abs(lat) > 60)
fig, ax, hm = heatmap(ncsst_polar[:,:,1])
ax.title[] = "regions north of 60°N and south of 60°S ($(ncsst_polar[:time][1]))"
Colorbar(fig[1,2],hm)
Makie.Colorbar()
Extract the time series of the closest point to a given point
ncsst_timeseries = @select(ncsst,lon ≈ 330 && lat ≈ 38)
lon = ncsst_timeseries["lon"][:]
lat = ncsst_timeseries["lat"][:]
nicetimeseries(ncsst_timeseries, title = "SST at $(lon) °E and $(lat) °N")
Extract the time series of the closest point with a specified tolerance If there is no data nearby an empty array is returned
ncsst_outside = @select(ncsst,lon ≈ 330 && lat ≈ 93 ± 0.1)
isempty(ncsst_outside)
size(ncsst_outside)
(0, 365)
Rather than selecting SST based on coordinates, we can also do the reverse: select the longitude based on SST. Find the longitudes where the SST exceeds 30 at the equator. sst
can contain missing values (e.g. on land). coalesce
is necessary here to replace the missing values by false in the boolean expression.
ds_equator = @select(ds,lat ≈ 0 && time ≈ DateTime(2023,1,1))
lon_equator = @select(ds_equator,coalesce(sst > 30,false))["lon"]
lon (54)
Datatype: Float32
Dimensions: lon
Attributes:
long_name = Longitude
standard_name = longitude
units = degrees_east
actual_range = Float32[0.125, 359.875]
axis = X
The first 3 longitude where SST exceeds 30°C at the equator for 2023-01-01
lon_equator[1:3]
3-element Vector{Float32}:
124.125
124.375
124.625
Grouping and reducing
With the function groupby
and macro @groupby
we can group the variable by a given criteria. For each group, we can then apply a reducting function (mean
, sum
, std
, ...).
For example group the SST by month and average per month:
sst_mean = mean(@groupby(ncsst,Dates.Month(time)))
sst (1440 × 720 × 12)
Datatype: Union{Missing, Float32}
Dimensions: lon × lat × time
Attributes:
long_name = Daily Sea Surface Temperature
units = degC
valid_range = Float32[-3.0, 45.0]
missing_value = -9.96921e36
precision = 2.0
dataset = NOAA High-resolution Blended Analysis
var_desc = Sea Surface Temperature
level_desc = Surface
statistic = Mean
parent_stat = Individual Observations
actual_range = Float32[-1.8, 36.94]
Instead of using the macro one can also use the function groupby
:
sst_mean = mean(groupby(ncsst,:time => Dates.Month))
nicemaps(sst_mean, timeindex = 1, title = "Mean January SST")
Use custom julia function for grouping: The Atlantic hurricane season is the period in a year, from June 1 through November 30, when tropical or subtropical cyclones are most likely to form in the North Atlantic Ocean. The data will be grouped into two cases: false
(time outside the Atlantic hurricane season) and true
(time is within the Atlantic hurricane season). We plot only the 2nd group.
sst_na = @select(ncsst,300 <= lon <= 360 && 0 <= lat <= 46)
sst_hs = mean(@groupby(sst_na,DateTime(year(time),6,1) <= time <= DateTime(year(time),11,30)))
nicemaps(sst_hs,timeindex = 2, title = "mean SST during the Atlantic hurricane season")
Use a user defined function is_atlantic_hurricane_season
without the @groupby
macro:
is_atlantic_hurricane_season(time) =
DateTime(year(time),6,1) <= time <= DateTime(year(time),11,30)
sst_std_hs = std(groupby(sst_na,:time => is_atlantic_hurricane_season));
nicemaps(sst_std_hs,timeindex = 2, title = "std SST during the Atlantic hurricane season")
Use a custom aggregation function to estimate the probability that the temperature exceeds 26°C during the Atlantic hurricane season
using CommonDataModel: GroupedVariable
prob_hot(x; dims=:) = mean(coalesce.(x .> 26,false); dims=dims)
prob_hot(gv::GroupedVariable) = reduce(prob_hot,gv)
is_atlantic_hurricane_season(time) =
DateTime(year(time),6,1) <= time <= DateTime(year(time),11,30)
sst_my_mean = prob_hot(groupby(ncsst,:time => is_atlantic_hurricane_season));
nicemaps(sst_my_mean,timeindex=2, title = "probability(temp. > 26°C) during Atlantic hurricane season")