Add support for downloading ERA5 pressure levels data product

examples/ERA5_bounded_pressure_level_data.jl

@@ -0,0 +1,121 @@
+# Download 3-D snapshots of the atmospheric state from ERA5, e.g., to be used with FieldTimeSeries
+#
+# ## Install dependencies
+#
+# ```julia
+# using Pkg
+# pkg"add NumericalEarth CDSAPI"
+# ```
+#
+# You also need CDS API credentials in `~/.cdsapirc`.
+# See <https://cds.climate.copernicus.eu/how-to-api> for setup instructions.
+#
+# See <https://cds.climate.copernicus.eu/datasets/reanalysis-era5-pressure-levels?tab=download> for
+# more information about pressure-level data
+
+using NumericalEarth
+using NumericalEarth.DataWrangling: Metadata, BoundingBox, download_dataset
+using NumericalEarth.DataWrangling.ERA5
+using CDSAPI
+using Dates
+using Oceananigans
+
+using Statistics
+using CairoMakie
+
+selected_levels = filter(≥(250hPa), ERA5_all_pressure_levels) # select all levels below 250 hPa
+dataset = ERA5HourlyPressureLevels(selected_levels)
+
+dates = DateTime(2004, 12, 16):Hour(1):DateTime(2005, 01, 09)  # van Zanten et al 2011 study period
+#dates = DateTime(2004, 12, 16):Hour(1):DateTime(2004, 12, 23)  # shorter time range for demo
+
+# Rauber et al 2007, Fig 1: focus region
+region = BoundingBox(latitude=(17, 18.5), longitude=(-62.5, -61))
+
+# OPTIONAL: download all variables at once (fewer CDS API requests)
+variables = [:geopotential, # for calculating zlevels as the mean geopotential height
+             :eastward_velocity,
+             :northward_velocity,
+             :temperature,
+             :specific_humidity]
+download_dataset(variables, dataset, dates; region)
+
+# ## Create time series
+
+T_meta = Metadata(:temperature;        dataset, dates, region)
+q_meta = Metadata(:specific_humidity;  dataset, dates, region)
+u_meta = Metadata(:eastward_velocity;  dataset, dates, region)
+v_meta = Metadata(:northward_velocity; dataset, dates, region)
+
+T_series = FieldTimeSeries(T_meta; time_indices_in_memory=2)
+q_series = FieldTimeSeries(q_meta; time_indices_in_memory=2)
+u_series = FieldTimeSeries(u_meta; time_indices_in_memory=2)
+v_series = FieldTimeSeries(v_meta; time_indices_in_memory=2)
+
+# ## Vertical profiles (mean ± spread)
+
+# Height coordinate (m) from the grid
+_, _, z = nodes(T_series[1])
+Nz = length(z)
+Nt = length(T_series.times)
+
+# Pressure at each level (hPa), ordered bottom-to-top (k=1 = highest pressure)
+p_levels = sort(selected_levels, rev=true) ./ hPa  # in hPa
+
+# Horizontal-mean profile for each snapshot → (Nz, Nt) arrays
+function horizontal_mean_profiles(series)
+    profiles = zeros(Nz, Nt)
+    for n in 1:Nt
+        data = interior(series[n], :, :, :)
+        profiles[:, n] = dropdims(mean(data, dims=(1, 2)), dims=(1, 2))
+    end
+    return profiles
+end
+
+T_profiles = horizontal_mean_profiles(T_series)
+q_profiles = horizontal_mean_profiles(q_series)
+u_profiles = horizontal_mean_profiles(u_series)
+v_profiles = horizontal_mean_profiles(v_series)
+
+# Convert T → potential temperature: θ = T (p₀/p)^(R/cₚ)
+Rₐ_over_cₚ = 0.286
+θ_profiles = T_profiles .* (1000 ./ p_levels) .^ Rₐ_over_cₚ
+
+# Convert specific humidity from kg/kg → g/kg
+q_profiles .*= 1000
+
+# ## Plot
+
+fig = Figure(size=(900, 500), fontsize=12)
+
+ax_θ = Axis(fig[1, 1], xlabel="θ [K]",       ylabel="Height [m]")
+ax_q = Axis(fig[1, 2], xlabel="qᵥ [g kg⁻¹]", ylabel="Height [m]")
+ax_u = Axis(fig[1, 3], xlabel="u [m s⁻¹]",   ylabel="Height [m]")
+ax_v = Axis(fig[1, 4], xlabel="v [m s⁻¹]",   ylabel="Height [m]")
+
+for (ax, profiles) in [(ax_θ, θ_profiles),
+                       (ax_q, q_profiles),
+                       (ax_u, u_profiles),
+                       (ax_v, v_profiles)]
+    μ = vec(mean(profiles, dims=2))
+    #lo = vec(minimum(profiles, dims=2))
+    #hi = vec(maximum(profiles, dims=2))
+    lo = [quantile(r, 0.25) for r in eachrow(profiles)]
+    hi = [quantile(r, 0.75) for r in eachrow(profiles)]
+    band!(ax, z, lo, hi; direction=:y, color=(:gray, 0.4))
+    lines!(ax, μ, z; color=:black, linewidth=2)
+end
+
+xlims!(ax_θ, 293, 317)
+xlims!(ax_q, 0, 20.4)
+xlims!(ax_u, -10, 0)
+xlims!(ax_v, -10, 0)
+
+linkyaxes!(ax_θ, ax_q, ax_u, ax_v)
+ylims!(ax_θ, 0, 4000)
+
+hideydecorations!(ax_q)
+hideydecorations!(ax_u)
+hideydecorations!(ax_v)
+
+save("ERA5_pressure_level_profiles.png", fig; px_per_unit=4)

examples/ERA5_bounded_surface_data.jl

@@ -0,0 +1,92 @@
+# Download 3-D snapshots of the atmospheric state from ERA5, e.g., to be used with FieldTimeSeries
+#
+# ## Install dependencies
+#
+# ```julia
+# using Pkg
+# pkg"add NumericalEarth CDSAPI"
+# ```
+#
+# You also need CDS API credentials in `~/.cdsapirc`.
+# See <https://cds.climate.copernicus.eu/how-to-api> for setup instructions.
+#
+# See <https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels?tab=download> for
+# more information about single-level data
+
+using NumericalEarth.DataWrangling: Metadata, BoundingBox, download_dataset
+using NumericalEarth.DataWrangling.ERA5
+using CDSAPI
+using Dates
+using Oceananigans
+
+using Statistics: mean
+using CairoMakie
+#precip_cmap = :rain
+precip_cmap = cgrad([:indigo, :darkblue, :blue, :deepskyblue, :cyan,
+                     :palegreen, :green, :yellow, :orange, :red, :pink])
+
+
+dataset = ERA5HourlySingleLevel()
+
+#dates = DateTime(2004, 12, 16):Hour(1):DateTime(2005, 01, 09)  # van Zanten et al 2011 study period
+dates = DateTime(2004, 12, 16):Hour(1):DateTime(2004, 12, 23)  # shorter time range for demo
+
+# Rauber et al 2007, Fig 1: precip map
+region = BoundingBox(latitude=(-25, 35), longitude=(-110, 30))
+
+# OPTIONAL: download all variables at once (fewer CDS API requests)
+# variables = [:total_precipitation,
+#              :sea_surface_temperature,
+#              :temperature, # at 2 m
+#              :dewpoint_temperature, # at 2 m
+#              :surface_pressure]
+# download_dataset(variables, dataset, dates; region)
+
+# ## Create time series
+
+precip_meta = Metadata(:total_precipitation; dataset, dates, region)
+precip_series = FieldTimeSeries(precip_meta)
+
+# ## Plotting
+
+Nt = length(precip_series.times)
+λ, φ, _ = nodes(precip_series[1])
+
+# ERA5 total_precipitation is in metres per hour; convert to mm/day
+to_mm_day = 1000 * 24
+
+# ### Time-averaged precipitation map
+
+precip_avg = mean(interior(precip_series[n], :, :, 1) for n in 1:Nt) .* to_mm_day
+
+fig1 = Figure(size=(900, 400))
+ax1 = Axis(fig1[1, 1],
+           title = "Mean precipitation (2004-12-16 to 2005-01-08)",
+           xlabel = "Longitude (°)",
+           ylabel = "Latitude (°)")
+ax1.xticks = -90:30:30
+hm = heatmap!(ax1, λ, φ, precip_avg; colormap=precip_cmap, colorrange=(0, 12))
+Colorbar(fig1[1, 2], hm, label="Precipitation (mm/day)")
+save("ERA5_mean_precipitation.png", fig1)
+
+# ### Precipitation animation
+
+fig2 = Figure(size=(900, 400))
+
+n = Observable(1)
+precip_n = @lift interior(precip_series[$n], :, :, 1) .* to_mm_day
+anim_title = @lift "Precipitation (mm/day), " * string(first(dates) + Second(round(Int, precip_series.times[$n])))
+
+ax2 = Axis(fig2[1, 1],
+           title = anim_title,
+           xlabel = "Longitude (°)",
+           ylabel = "Latitude (°)")
+ax2.xticks = -90:30:30
+
+hm2 = heatmap!(ax2, λ, φ, precip_n; colormap=precip_cmap, colorrange=(0, 12))
+Colorbar(fig2[1, 2], hm2, label="Precipitation (mm/day)")
+
+record(fig2, "ERA5_precipitation.mp4", 1:Nt; framerate=12) do nn
+    n[] = nn
+end
+

examples/ERA5_winds_and_stokes_drift.jl

@@ -16,7 +16,7 @@
 
 using NumericalEarth
 using NumericalEarth.DataWrangling: Metadatum
-using NumericalEarth.DataWrangling.ERA5: ERA5Hourly
+using NumericalEarth.DataWrangling.ERA5: ERA5HourlySingleLevel
 using CDSAPI
 
 using Oceananigans
@@ -29,7 +29,7 @@ using Dates
 # while ocean wave variables (Stokes drift) live on a 0.5° grid (720×361).
 # We define metadata for each variable at a single date.
 
-dataset = ERA5Hourly()
+dataset = ERA5HourlySingleLevel()
 date = DateTime(2020, 1, 15, 12) # January 15, 2020 at 12:00 UTC
 
 u_stokes_meta = Metadatum(:eastward_stokes_drift;  dataset, date)
@@ -42,10 +42,11 @@ v_wind_meta   = Metadatum(:northward_velocity;     dataset, date)
 # We build a single `LatitudeLongitudeGrid` and use `set!` to download
 # and interpolate all four variables onto it.
 
-grid = LatitudeLongitudeGrid(size = (1440, 721, 1),
+grid = LatitudeLongitudeGrid(size = (1440, 720, 1),
                              longitude = (0, 360),
                              latitude = (-90, 90),
                              z = (0, 1))
+@show grid
 
 u_stokes = CenterField(grid)
 v_stokes = CenterField(grid)