Using the xarray structured data interface within Salvus

Originally, Salvus entirely relied on obspy.Stream objects to internally represent simulated and observed data. While obspy.Stream’s are a great representation type for seismological data and come with convenient processing and filtering functions, there are a few limitations to the original data interface used within Salvus.

• Data could only be retrieved per station. If an event contained multiple receivers, the data for each receiver had to be retrieved separately in a loop, which could become slow in cases where a high number of receivers was used.
• All data processing, windowing, and the computation of misfits was performed trace-by-trace. Besides the already-mentioned performance aspects, this also meant that more sophisticated time-and-space-dependent processing and misfit computations (e.g. in the fk-domain) were not supported.
• obspy.Stream objects are mainly used within the seismological community. For non-seismology users, learning how to work with these objects would provide an unnecessary extra challenge.

For these reasons, with salvus release 2025.1.0, we have changed the way that data is represented internally within Salvus. Salvus now uses xarray data structures with a specific format to represent data internally. In this tutorial, you will learn what this format looks like and how you can start using it within Salvus. If you prefer to continue to use the old interface based on obspy.Stream objects, do not worry, all the old functionality is still supported!

This tutorial is structured in the following way:

1. Introduction to xarray (if you are already familiar with xarray, you can skip this step).
2. Retrieval of xarray data from a Salvus project.
3. Custom data processing, data selection, and misfit functions using Salvus’ xarray interface.

1. Introduction to xarray

Xarray provides a convenient interface to represent structured data. In a nutshell, an xarray data structure is simply a multidimensional array with a label for each dimension. It can optionally store additional information on the coordinates of each data point. Additionally, xarray comes with some convenient features like data selection, plotting, and sorting methods.

To familiarize ourselves with the xarray interface, let us generate a simple 2D data set analogous to a seismic data set. In this simple example, we model a space- and time-dependent seismic wavefield with the following simple monochromatic plane-wave solution to the wave equation

$u(x, t) = \sin(\omega t - k x)$,

with $x$ being the space coordinate, $t$ the time coordinate, $\omega$ the angular frequency, and $k$ the wavenumber, which can be expressed in terms of the velocity $c$, as $k=\omega/c$.

Let us use this function to generate a discrete seismic data set for a monochromatic plane wavefield at 10 Hz where we sample the wavefield along $x$ with a linear array of 25 sensors spaced at 10 meters. The time axis is resolved with 20 samples per period.

Note that we obtained our data as a 2D numpy array, which does not contain any information on the time and the space coordinates. Here, xarray can help us to additionally store this information in a single data structure. Let us therefore wrap our data into an xarray.DataArray structure. Upon creation of the xarray object can provide the following information:

• data: The numpy array containing our data.
• dims: A list with the labels for each dimension of our dataset. Our data is 2D, so we add two labels, one for the time and one for the space dimension.
• coords: Here, we can optionally provide the coordinates for each dimension. We pass the time and space axis as a dictionary. The keys of the dictionary have to agree with the dimensions the we set in the previous row.
• name: An optional name for our data.

Once we have stored our data in an xarray.DataArray, we can use the plot() method for a quick visualization. Note that the time and space coordinates are now represented correctly, as well as adding an appropriate color scale automatically.

We can easily retrieve the underlying data as numpy arrays in the following way using the data attribute:

Coordinates do not necessarily need to be numerical values. For example, instead of the x-coordinate, we might want to store something like a receiver name that uniquely identifies the receiver that recorded the data.

To select data of a specific receiver, we can use the sel() method on the "receiver_id" dimension:

For index based selection, we can use the isel() method:

Finally, it might be interesting to know that one can construct Datasets from multiple DataArrays. These objects can be useful when dealing with multiple quantities together, such as medium models or multiple recording fields.

Multi-index dimensions

In certain cases, our data will have more dimensions than just the space and time dimension. For example, elastic seismic wavefields are vector fields of which we might record multiple components. Our data will thus have an additional dimension corresponding to the component of the wavefield. For various reasons, data within Salvus are always represented as 2D xarray data structures. In order to keep track of information on additional dimensions such as the component, we make use a so-called multi-index dimension that can store information on both the "receiver_id" and the "component". Let us see how this works in practice.

First, we create a 3D data array, where the first dimension corresponds to the "receiver_id", the second dimension to the "component", and the third dimension to "time".

Now, we set up a 3D xarray.DataArray from that.

We now convert this 3D array to a 2D array by combining the "receiver_id" and "component" dimensions into a multi-indexed "trace" dimension. To do so, we use the stack() method. This is exactly how data is represented internally within Salvus: any dimension that is not time is stacked into trace.

We can still select data of one of the original dimensions, like a single component, the same way:

Or, we can directly make a selection on the new combined "trace" dimension for a specific receiver and component pair like this.

This won’t work, as it is not the correct order of coordinates. da_multi_indexed.sel(trace=(“Y”, “R11”)).plot()

To go back to a 3D representation of the data, we can use the unstack method, which gives us our original 3D array.

2. Retrieval of xarray data structures within Salvus

Now that we know some of the basics of xarray data structures, let us have a closer look at how we can use them within Salvus. In the following example, we use Salvus to simulate some simple 2D seismic line data as it would be recorded by a line of geophones at the Earth’s surface. The following cell sets up a simple elastic 2D simulation in a homogeneous domain with the following elastic properties:

• P-wave velocity = 3000 m/s
• S-wave velocity = 1500 m/s
• Density = 2000 kg/m$^3$

We set up a single vertically-directed source and a line of receivers spaced at 10 meters. Please refer to the tutorial “My first Salvus simulation” for more details on how to set up a Salvus simulation.

Now we run the simulations.

So, how do we now access the data? First, let us revisit how we can retrieve data of a single receiver as an obspy stream object. To do so, we first retrieve an EventData object from which we can then get the data using the get_waveform_data() method. For example, to retrieve data for the first receiver with id XX.AA_1., we would do the following.

If we want to avoid retrieving the data for each receiver separately by looping over all receiver names, we can directly access all data of this event as an xarray.DataArray using the function get_waveform_data_xarray() that was introduced with Salvus release 2025.1.0.

Looking at the output, we observe the following structure of the data.

The data is wrapped in a 2D array with the name corresponding to the field of the data, in this case "displacement". The array has the following two dimensions, always in this order:

• "trace": A multi-index dimension that combines the dimensions "receiver_id" and "component". For details on multi-indices, see section 1 of this tutorial.
• "time: The time dimension.

The data contains additional attributes. The most important ones are the following:

• "start_time_in_seconds": The time of the first sample.
• "sampling_rate_in_hertz": The sampling rate of the data.
• "npts": The number of time samples.
• "location": The coordinates of each receiver location.

We can now in turn select the two components "X" and "Y" and plot them.

Again, we can use the unstack() method to obtain a 3D representation of the data with dimensions ("receiver_id", "component", "time").

3. Custom processing, selection, and misfits using xarrays

You can use SalvusProject to apply custom processing, data selection, and misfit functions to you data. To do so, you need to format your custom functions in a specific way. As of Salvus release 2025.1.0, you can now define custom functions that operate on the new xarray data structures directly and add them to a project. Trace-by-trace based custom functions operating on obspy.Stream objects are still supported as well. However, you should observe a significant performance improvement when using the xarray interface for custom processing functions, especially when working with data with a large number of receivers.

In the following, we briefly describe the format of the new custom functions and give some examples on how you can use these functions within a Salvus project.

Data processing / filtering

First, let us consider the case where we want to apply some filtering to our data. Here, we consider a simple bandpass filter that filters the data to a frequency band beteween 10 and 20 Hz. This function will directly operate on an xarray.DataArray object.

If we want Salvus to recognize our custom filtering function, we need to specify it in a specific way. The function needs to take exactly two mandatory input arguments:

• data_array: The Salvus xarray.DataArray object on which the filtering will be performed.
• event: The Salvus Event object of the Event that this data corresponds to.

Let us apply this function to our simulated data to see if this works as intended.

This looks as expected.

We can now directly add this function to our project so that Salvus can use it internally. This is achieved via a ProcessingConfiguration.

The project is now aware of this function and we can use it, for example, when defining a new misfit configuration.

To retrieve filtered data from the project, we can access it in the following way, adding the function “bandpass” with a : to the data that is being queried:

Data selection / windowing

In a similar fashion, we can define custom data selection and windowing functions using the xarray interface. Here, we want to define a data selection function that performs a couple of operations on our input data:

First, the data selection function should only select a subset of the available receivers. This is important if we want to disregard certain receivers in an inversion (e.g., receivers that record data that is too noisy).

Second, the function should only select a single component. In this case, we want to select only the vertical component.

Last, the data selection function should apply a time window to the data.

Again, our data selection function will need to accept two mandatory input arguments:

• data_array: The Salvus xarray.DataArray object on which the data selection will be performed.
• event: The Salvus Event object of the event that this data corresponds to.

The output of the function should then again be a DataArray, with the same dimensions as the DataArray on which the selection was performed. However, please note that in this case, not all coordinates of the original DataArray need to be present: the omission of e.g. a receiver_id or component implicitly sets its weights to zero. The data of the weights array itself should correspond to the temporal weighting function (window) that will be applied to the data.

In the following, we define such a data selection function.

We can now create a custom DataSelectionConfiguration and add it to our project.

This data selection configuration can now be used when setting up a MisfitConfiguration. To check whether the data selection works as intended, we retrieve data from the project and specify that we want to use our new data selection configuration.

Note that the resulting DataArray has the temporal weights applied, but any station or component that wasn’t present in the returned object from select_data is not present at all in the retrieved data.

If we now plot the output, you will see that our data selection works as intended.

Custom misfit functions

In earlier versions of Salvus, misfits were always computed on a trace-by-trace basis by looping over all receivers and components. Now, misfits are directly computed on the new xarray structures, which allows us to vectorize the misfit computations. If you want to define a custom misfit function, you can now take advantage of this. Additionally, it is now possible to define time-and-space dependent misfit functions.

Here, we simply review the standard L2 misfit to understand how custom misfit functions can be implemented within Salvus.

A custom misfit function needs to take three mandatory input arguments:

• data_array_synthetic: The synthetic waveforms as an xarray.DataArray object.
• data_array_observed: The observed waveforms as an xarray.DataArray object.
• event: The event object.

The output must be an xarray.Dataset object. A xarray.Dataset object combines multiple xarray.DataArray’s. The output Dataset will contain two DataArray’s, one for the adjoint sources and one for the per-trace misfits.

The following code snippets demonstrates the computation of the misfits and adjoint sources for a conventional L2 misfit.

We can now set up a new misfit configuration using all of our custom functions (processing, data_selection, and custom misfit) and add it to our project in the following way.

Let us verify that everything was properly added to our project: