using Earth2Studio
using Dates
using DimensionalData
using CairoMakie
using PythonCall: pyimport, pyconvert, pylist, pydictJulia packages
earth2studio returns xarray objects from data sources and IO backends. DimensionalData.jl converts xarray.DataArray and xarray.Dataset to DimArray and DimStack without copying the data, and its Makie extension plots the result with labelled axes. The PythonCall integration activates automatically when both packages are loaded.
xarray to DimensionalData
Fetch 2 m temperature and 10 m zonal wind at five times from ARCO ERA5:
ds = Earth2Studio.data.ARCO(verbose=false)
t0 = DateTime(2023, 6, 15)
da = ds(pylist(t0 .+ Hour.(0:6:24)), pylist(["t2m", "u10m"]))
da.dimsPython: ('time', 'variable', 'lat', 'lon')pyconvert(DimArray, da) does the conversion, with one preparation step. DimensionalData reads each coordinate as an attribute of the DataArray, and DataArray.variable is an xarray property that returns the array’s underlying Variable rather than the coordinate named variable. Rename that dimension before converting:
A = pyconvert(DimArray, da.rename(variable="field"))┌ 1440×721×2×5 DimArray{Float64, 4} ┐ ├───────────────────────────────────┴──────────────────────────────────── dims ┐ ↓ lon Sampled{Float64} [0.0, …, 359.75] ForwardOrdered Irregular Points, → lat Sampled{Float64} [90.0, …, -90.0] ReverseOrdered Irregular Points, ↗ field Categorical{PythonCall.Utils.StaticString{UInt32, 4}} ["t2m", "u10m"] ForwardOrdered, ⬔ time Sampled{PythonCall.NumpyDates.InlineDateTime64{PythonCall.NumpyDates.MICROSECONDS}} [PythonCall.NumpyDates.InlineDateTime64{PythonCall.NumpyDates.MICROSECONDS}("2023-06-15T00:00:00"), …, PythonCall.NumpyDates.InlineDateTime64{PythonCall.NumpyDates.MICROSECONDS}("2023-06-16T00:00:00")] ForwardOrdered Irregular Points └──────────────────────────────────────────────────────────────────────────────┘ [:, :, 1, 1] ↓ → 90.0 89.75 89.5 … -89.25 -89.5 -89.75 -90.0 0.0 273.981 273.941 273.922 224.077 224.429 224.98 224.723 0.25 273.981 273.941 273.921 224.079 224.43 224.985 224.723 0.5 273.981 273.941 273.921 224.08 224.434 224.989 224.723 ⋮ ⋱ ⋮ 359.25 273.981 273.941 273.922 224.079 224.424 224.977 224.723 359.5 273.981 273.941 273.922 224.079 224.427 224.979 224.723 359.75 273.981 273.941 273.922 … 224.077 224.429 224.98 224.723
Two things changed. The dimensions are reversed, (lon, lat, field, time) in place of (time, variable, lat, lon), because Python stores arrays row-major and Julia column-major. Reversing keeps the fast axis the same in both languages, so no data moves. And the data is shared: parent(A) is a PyArray viewing the numpy buffer. Pass copy=true to get a plain Array.
typeof(parent(A))PythonCall.PyArray{Float64, 4, true, true, Float64}The time lookup holds numpy datetime64 values. Replace it with DateTime so that selectors take Julia dates and Makie can label the axis:
A = set(A, :time => pyconvert(Vector{DateTime}, da.time.values))
lookup(A, :time)Sampled{Dates.DateTime} ForwardOrdered Irregular DimensionalData.Dimensions.Lookups.Points
wrapping: 5-element Vector{DateTime}:
2023-06-15T00:00:00
2023-06-15T06:00:00
2023-06-15T12:00:00
2023-06-15T18:00:00
2023-06-16T00:00:00
Selectors index by coordinate value. At matches exactly and Near picks the closest point:
raleigh = A[field=At("t2m"), lat=Near(35.8), lon=Near(281.3)]┌ 5-element DimArray{Float64, 1} ┐ ├────────────────────────────────┴─────────────────────────────────────── dims ┐ ↓ time Sampled{Dates.DateTime} [Dates.DateTime("2023-06-15T00:00:00"), …, Dates.DateTime("2023-06-16T00:00:00")] ForwardOrdered Irregular Points └──────────────────────────────────────────────────────────────────────────────┘ 2023-06-15T00:00:00 300.231 2023-06-15T06:00:00 294.853 2023-06-15T12:00:00 294.244 2023-06-15T18:00:00 302.782 2023-06-16T00:00:00 300.388
Plotting with Makie
DimensionalData’s Makie extension reads axis labels, ticks, and the title from the dimensions, so a DimArray plots without extracting coordinate vectors. The title lists the coordinates that were selected away.
fig = Figure(size=(900, 450))
heatmap(fig[1, 1], A[field=At("t2m"), time=At(t0)] .- 273.15; colormap=:thermal)
fig
lines(raleigh .- 273.15; axislegend=false)
Datasets to DimStack
XarrayBackend stores a forecast as an xarray.Dataset with one data variable per field. pyconvert(DimStack, ...) converts each variable to a layer that shares the same dimensions. The setup below is the Persistence forecast from the Forecasting tutorial.
np = pyimport("numpy")
grid = pydict("lat" => np.linspace(90, -90, 721), "lon" => np.linspace(0, 359.75, 1440))
model = Earth2Studio.models.px.Persistence("t2m", grid)
io = Earth2Studio.io.XarrayBackend()
Earth2Studio.run.deterministic([t0], 4, model, ds, io; verbose=false)
S = pyconvert(DimStack, io.root)┌ 1440×721×5×1 DimStack ┐ ├───────────────────────┴──────────────────────────────────────────────── dims ┐ ↓ lon Sampled{Float64} [0.0, …, 359.75] ForwardOrdered Irregular Points, → lat Sampled{Float64} [90.0, …, -90.0] ReverseOrdered Irregular Points, ↗ lead_time Sampled{PythonCall.NumpyDates.InlineTimeDelta64{PythonCall.NumpyDates.SECONDS}} [PythonCall.NumpyDates.InlineTimeDelta64{PythonCall.NumpyDates.SECONDS}(0), …, PythonCall.NumpyDates.InlineTimeDelta64{PythonCall.NumpyDates.SECONDS}(86400)] ForwardOrdered Irregular Points, ⬔ time Sampled{PythonCall.NumpyDates.InlineDateTime64{PythonCall.NumpyDates.NANOSECONDS}} [PythonCall.NumpyDates.InlineDateTime64{PythonCall.NumpyDates.NANOSECONDS}("2023-06-15T00:00:00")] ForwardOrdered Irregular Points ├────────────────────────────────────────────────────────────────────── layers ┤ :t2m eltype: Float64 dims: lon, lat, lead_time, time size: 1440×721×5×1 └──────────────────────────────────────────────────────────────────────────────┘
The lead_time lookup holds numpy timedelta64 values, which convert to Hour the same way time converts to DateTime. Layers are accessed by name and support the same selectors:
S = set(S, :time => pyconvert(Vector{DateTime}, io.root.time.values),
:lead_time => pyconvert(Vector{Hour}, io.root.lead_time.values))
S[:t2m][lead_time=At(Hour(12)), lat=Near(35.8), lon=Near(281.3)]┌ 1-element DimArray{Float64, 1} t2m ┐ ├────────────────────────────────────┴─────────────────────────────────── dims ┐ ↓ time Sampled{Dates.DateTime} [Dates.DateTime("2023-06-15T00:00:00")] ForwardOrdered Irregular Points └──────────────────────────────────────────────────────────────────────────────┘ 2023-06-15T00:00:00 300.231
Other packages
Rasters.jl and YAXArrays.jl build on DimensionalData, so a converted DimArray is their input type once its lon and lat dimensions are renamed to X and Y. For a Python-free path, NetCDF4Backend writes a file that NCDatasets.jl and Rasters.jl open directly.