workingWith3Ddatasets

Loading a netcdf dataset and working with xarray¶

In this notebook, we explore a 3D NetCDF dataset of seismic anisotropy derived from the Eastern Mediterranean sector of the Collaborative Seismic Earth Model (CSEM). This high-resolution tomographic model offers valuable insights into the anisotropic structure of the Earth’s crust and mantle in a tectonically complex region.¶

The dataset is publicly available through the IRIS Earth Model Collaboration at: 🔗 https://ds.iris.edu/ds/products/emc-csem_eastmed/

We will demonstrate how to:

Load and inspect the 3D NetCDF data

Visualize seismic anisotropy parameters in 2D slices and 3D volumes

Extract physical insights using Python-based geophysical toolkits

This work supports broader efforts to understand lithospheric deformation and mantle flow in one of the most geodynamically active regions on Earth.

Library needed¶

xarray: xarray is a powerful Python library designed for working with labeled multi-dimensional arrays. It builds on NumPy and pandas, adding support for dimensions (e.g., latitude, longitude, depth, time) and coordinates, making it particularly suited for working with NetCDF files and other scientific datasets.

In [1]:

import xarray as xr
import numpy as np

To open a dataset and store it in a variable, we use open_dataset() method of xarray

In [2]:

ds = xr.open_dataset("/Volumes/work/work/data/csem-eastmed-2019.12.01.nc")

This is how ds looks:

Note: The dimension of the dataset is 7716681 which is how all data variables will be written

To list the data variables available, we use ds.data_vars

In [4]:

ds.data_vars

Out[4]:

Data variables:
    vsv      (depth, latitude, longitude) float32 4MB ...
    vsh      (depth, latitude, longitude) float32 4MB ...

Coordinates are the dimensions of the dataset, which can be accessed using ds.coords:

In [5]:

ds.coords

Out[5]:

Coordinates:
  * latitude   (latitude) float32 308B 30.0 30.2 30.4 30.6 ... 44.8 45.0 45.2
  * longitude  (longitude) float32 664B 12.0 12.2 12.4 12.6 ... 44.6 44.8 45.0
  * depth      (depth) float32 324B 0.0 5.0 10.0 15.0 ... 390.0 395.0 400.0

All other metadata related to this dataset is stored in attrs and can be accessed by ds.attrs

In [6]:

ds.attrs

Out[6]:

{'title': 'Collaborative Seismic Earth Model (CSEM) - extraction of the Eastern Mediterranean and Turkey',
 'id': 'csem_eastmed',
 'summary': 'CSEM extraction of the Eastern Mediterranean and Turkey',
 'keywords': 'Collaborative Seismic Earth Model, multi-scale full-waveform inversion',
 'Conventions': 'CF-1.0',
 'Metadata_Conventions': 'Unidata Dataset Discovery v1.0',
 'acknowledgement': 'provided by ETH SWP Group, see publication references',
 'history': '2020-10-08 13:11:04 UTC Converted to netCDF by GeoCSV_2_netCDF_3D.py V.2020.273 from csem-eastmed-2019.12.01.csv\n2020-10-08: IRIS DMC updated reference_pid; 2020-03-02 19:34:09 UTC Converted to GeoCSV by netCDF_2_GeoCSV_3D.py ,V.2020.006 from csem-eastmed-2019.12.01.nc; created Fri Nov 22 16:42:25 2019; 2020-02-28 IRIS DMC, updated metadata to organize reference, author, repository and also add PID',
 'comment': '',
 'geospatial_lat_min': 30.0,
 'geospatial_lat_max': 45.0,
 'geospatial_lat_units': 'degrees_north',
 'geospatial_lat_resolution': 0.2,
 'geospatial_lon_min': 12.0,
 'geospatial_lon_max': 45.0,
 'geospatial_lon_units': 'degrees',
 'geospatial_lon_resolution': 0.2,
 'geospatial_vertical_min': 0.0,
 'geospatial_vertical_max': 400.0,
 'geospatial_vertical_units': 'km',
 'geospatial_vertical_positive': 'down',
 'netcdf_file': 'csem-eastmed-2019.12.01.nc',
 'reference': 'Fichtner, Cubuk-Sabuncu, Blom, and Gokhberg (2019)',
 'reference_pid': 'doi:10.5194/se-11-669-2020; doi:10.1016/j.pepi.2017.06.014; doi:10.1093/gji/ggt118; doi:10.1029/2018GL077338',
 'author_name': 'Andreas Fichtner',
 'author_email': 'andreas.fichtner@erdw.ethz.ch',
 'author_institution': 'ETH Zurich, Zurich, Switzerland',
 'author_url': 'https://cos.ethz.ch/research/CSEM.html',
 'repository_name': 'EMC',
 'repository_institution': 'IRIS DMC',
 'repository_pid': 'doi:10.17611/dp/emccsemeastmed20191201',
 'model': 'CSEM_Eastmed'}

In [7]:

ds.data_vars.keys

Out[7]:

<bound method Mapping.keys of Data variables:
    vsv      (depth, latitude, longitude) float32 4MB ...
    vsh      (depth, latitude, longitude) float32 4MB ...>

In [8]:

ds.data_vars

Out[8]:

Data variables:
    vsv      (depth, latitude, longitude) float32 4MB ...
    vsh      (depth, latitude, longitude) float32 4MB ...

The dataset contains two key velocity variables:¶

vsv: Vertically polarized S-wave velocity

vsh: Horizontally polarized S-wave velocity

To compute radial anisotropy—which quantifies the difference in shear wave speeds due to anisotropic structure—we use the following formula:

$$ \text{Radial Anisotropy (\%)} = 100 \times \frac{v_{sh} – v_{sv}}{\left(\frac{v_{sh} + v_{sv}}{2}\right)} $$

In [9]:

radial_anisotropy = 100 * (ds['vsh'] - ds['vsv']) / ((ds['vsh'] + ds['vsv']) / 2)

The upper code creates yet another xarray with dimension depth:81*latitude:77*longitude:66

To merge this into the earlier array we use this command

In [11]:

ds['radial_anisotropy'] = radial_anisotropy

If we want to add a metadata about radial_anisotropy, we can do this

In [13]:

ds['radial_anisotropy'].attrs['long_name'] = 'Radial Anisotropy'
ds['radial_anisotropy'].attrs['units'] = '%'
ds['radial_anisotropy'].attrs['description'] = '100 * (vsh - vsv) / ((vsh + vsv) / 2)'

Lets add another variable, absolute radial anisotropy:

$$ \text{Absolute Radial Anisotropy (\%)} = abs(100 \times \frac{v_{sh} – v_{sv}}{\left(\frac{v_{sh} + v_{sv}}{2}\right)}) $$

In [14]:

absolute_radial_anisotropy = abs(100 * (ds['vsh'] - ds['vsv']) / ((ds['vsh'] + ds['vsv']) / 2))

In [15]:

ds['absolute_radial_anisotropy'] = absolute_radial_anisotropy

To save this dataset with the new variable, we can use the to_netcdf method:

In [17]:

ds.to_netcdf("/Volumes/work/utpalsingh.in/csem-eastmed-2019.12.01-with-radial-anisotropy.nc")

Creating a new netcdf dataset¶

Lets create a new dataset and save to netcdf with just radial_anisotropy and absolute radial anisotropy

In [18]:

ds_aniso = xr.Dataset(
    data_vars={
        'radial_anisotropy': (['depth', 'latitude', 'longitude'], radial_anisotropy.data),
        'absolute_radial_anisotropy': (['depth', 'latitude', 'longitude'], absolute_radial_anisotropy.data)
    },
    coords={
        'depth': ds.depth,
        'latitude': ds.latitude,
        'longitude': ds.longitude
    },
    attrs={
        'title': 'Radial Anisotropy Dataset',
        'source': 'Derived from CSEM vsv and vsh',
        'description': 'Contains only the computed radial anisotropy in percentage (%)'
    }
)

In [20]:

ds_aniso.to_netcdf("/Volumes/work/utpalsingh.in/csem-eastmed-2019.12.01-radial-anisotropy.nc")

Lets plot radial anisotropy using pyvista

In [21]:

import pyvista as pv
import numpy as np

# Extract data
data = ds.radial_anisotropy.values  # Shape: (depth, lat, lon)

# Get the physical coordinates
depth = ds.depth.values
lat = ds.latitude.values
lon = ds.longitude.values

data = data.astype(np.float32)

# Get grid spacing
dz = np.diff(depth).mean()
dy = np.diff(lat).mean()
dx = np.diff(lon).mean()

# Build the grid dimensions (X fastest, Z slowest)
nx, ny, nz = len(lon), len(lat), len(depth)
spacing = (dx, dy, dz)

# Create a ImageData
grid = pv.ImageData()

grid.dimensions = (nx, ny, nz)  # Note: +1 not needed unless you have node-centered data
grid.origin = (lon[0], lat[0], depth[0])  # (x0, y0, z0)
grid.spacing = spacing  # (dx, dy, dz)

# Transpose data to match PyVista ordering (z, y, x)
data = np.transpose(data, (2, 1, 0))  # from (depth, lat, lon) to (lon, lat, depth)

# Flatten and assign data to grid
grid["radial_anisotropy"] = data.ravel(order="F")  # Fortran order matches grid memory layout

# Plot using volume rendering
p = pv.Plotter()
# p.add_volume(grid, scalars="radial_anisotropy", opacity="sigmoid", cmap="coolwarm", clim=(-12, 12))
p.add_volume(grid, scalars="radial_anisotropy", cmap="coolwarm", clim=(-12, 12))

p.set_scale(xscale=1.0, yscale=1.0, zscale=0.02)

p.add_axes()
p.show()

2025-07-29 22:51:45.469 (   0.929s) [           80054]        vtkRenderer.cxx:1161  WARN| vtkOpenGLRenderer (0x3040ef200): Resetting view-up since view plane normal is parallel

Widget(value='<iframe src="http://localhost:56394/index.html?ui=P_0x17cf9d580_0&reconnect=auto" class="pyvista…

This will generate a plot like this

To make plots using matplotlib, we do this¶

In [25]:

# Plot horizontal slice at depth 50km

import matplotlib.pyplot as plt
import numpy as np

target_depth = 50.0
depth_index = np.abs(ds.depth.values - target_depth).argmin()
print(depth_index)

# Extract the 2D slice
anisotropy_slice = ds.radial_anisotropy.isel(depth=depth_index)

# Plot
plt.figure(figsize=(10, 6))
im = plt.pcolormesh(
    ds.longitude, ds.latitude, anisotropy_slice,
    shading='auto', cmap='coolwarm', vmin=-12, vmax=12
)
plt.colorbar(im, label='Radial Anisotropy (%)')
plt.title(f'Radial Anisotropy at {ds.depth.values[depth_index]:.1f} km Depth')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.grid(True)
plt.tight_layout()
plt.show()

No description has been provided for this image

In [30]:

Out[30]:

array([ 50.,  45.,  40.,  35.,  30.,  25.,  20.,  15.,  10.,   5.,   0.,
         5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.,  50.,  55.,
        60.,  65.,  70.,  75.,  80.,  85.,  90.,  95., 100., 105., 110.,
       115., 120., 125., 130., 135., 140., 145., 150., 155., 160., 165.,
       170., 175., 180., 185., 190., 195., 200., 205., 210., 215., 220.,
       225., 230., 235., 240., 245., 250., 255., 260., 265., 270., 275.,
       280., 285., 290., 295., 300., 305., 310., 315., 320., 325., 330.,
       335., 340., 345., 350.], dtype=float32)

In [32]:

# Plot a Vertical Cross-Section at a Fixed Longitude

# Choose a longitude to slice
target_lon = 30.0
lon_index = np.abs(ds.longitude.values - target_lon).argmin()
print(lon_index)

# Extract 2D cross-section (depth vs latitude)
section = ds.radial_anisotropy.isel(longitude=lon_index)

# Plot
plt.figure(figsize=(10, 6))
im = plt.imshow(
    section.T,  # Transpose for correct orientation: (depth, lat)
    cmap='coolwarm', vmin=-12, vmax=12,
    aspect='auto',
    origin='upper',
    extent=[ds.latitude.min(), ds.latitude.max(), ds.depth.max(), ds.depth.min()]
)
plt.colorbar(im, label='Radial Anisotropy (%)')
plt.title(f'Radial Anisotropy Cross-Section at {ds.longitude.values[lon_index]:.1f}°E')
plt.xlabel('Latitude')
plt.ylabel('Depth (km)')
plt.grid(True)
plt.tight_layout()
plt.show()

In [ ]:

Working with 3D Geo-Datasets and plotting

Loading a netcdf dataset and working with xarray¶

Library needed¶

The dataset contains two key velocity variables:¶

Creating a new netcdf dataset¶

To make plots using matplotlib, we do this¶

Leave a Comment Cancel Reply