---
title: volume rendering intro for geodata in yt
tags: yt
---

# volume rendering intro for geodata in yt

this presentation: https://bit.ly/ytgeostuff
yt's overview on volume rendering: [link](https://yt-project.org/doc/visualizing/volume_rendering.html)
Chris's 2020 AGU poster: [citeable archive](https://www.essoar.org/doi/abs/10.1002/essoar.10506118.2), [direct repo link](https://github.com/earthcube2020/ec20_havlin_etal).

Other code/repos:

yt: https://yt-project.org/
ytgeotools: https://github.com/chrishavlin/ytgeotools
yt_idv: https://github.com/yt-project/yt_idv

## where to get seismic data

tomography http://ds.iris.edu/ds/products/emc-earthmodels/

## general workflow for seismic data

First step is to read 3d model data from native storage (e.g., netcdf file) into a yt in-memory dataset `yt.load_uniform_grid()` [docs link](https://yt-project.org/doc/reference/api/yt.loaders.html#yt.loaders.load_uniform_grid)

How to do that depends: making maps or 3D renderings?

### map making

load as an `internal_geographic` dataset:

ds = yt.load_uniform_grid(data,  
                          sizes, 
                          1.0, 
                          geometry=("internal_geographic", dims),
                          bbox=bbox)

where:

`data`: dictionary of field data arrays read from your netcdf
`sizes`: the shape of the field data
`bbox`: the bounding box
`dims`: the ordering of dimensions in the data arrays

See [this notebook](https://nbviewer.jupyter.org/github/chrishavlin/AGU2020/blob/main/notebooks/seismic.ipynb#Fixed-Depth-Maps) for a full example.

See [this notebook](https://nbviewer.jupyter.org/github/chrishavlin/yt_scratch/blob/master/notebooks/yt_obspy_raypath_sampling.ipynb) for an example of some analysis once data is loaded: sampling data along a seismic ray path.

Caveats:
* fixed-depth maps only (no cross-sections)
* still a little buggy (global models seem to work better than non-global models)
* no volume rendering

## geo-volume rendering with yt (background)

**data must be in cartesian coordinates**

Step 1: interpolate the seismic model to a (uniform) cartesian grid.

My approach:

1. convert lat, lon, depth arrays to geocentric cartesian coordinates, `x,y,z`
2. find bounding box in cartesian coordinates:

cart_bbox = [
             [min_x, max_x], 
             [min_y, max_y],
             [min_z, max_z],
            ]

3. create a uniform grid to cover bounding box:

x_i = np.linspace(min_x, max_x, n_x) 
y_i = np.linspace(min_y, max_y, n_y) 
z_i = np.linspace(min_z, max_z, n_z) 

(`i` for interpolation)

4. find the data values at all those points

data_on_cart_grid = new array
for xv, yv, zv in np.meshgrid(x_i, y_i, z_i):
    data_on_cart_grid[xvi, yvi, zvi] = sample_the_data(xv, yv, zv)

**How to `sample_the_data` ?**

Build a KDtree of actual model points

tree = KDtree(x, y, z)

for every point on the new cartesian grid, find the `N` closest points, use inverse-distance weighting (IDW) to average those points. Can set a max distance to exclude points too far away.

Could experiment with different methods here! just take the "nearest"? Different interopolation method entirely?

Sometimes a little buggy... depending on the resolution of the covering grid: weird artifacts or empty interpolated data.

**Once data is sampled**

Use `yt.load_uniform_grid`:

data = {"dvs": dvs_on_cartesian}
bbox = cart_bbox
sizes = dvs_on_cartesian.shape
ds = yt.load_uniform_grid(data, sizes, 1.0)

can now use volume rendering!

Seismic examples:
* [AGU2020](https://nbviewer.jupyter.org/github/chrishavlin/AGU2020/blob/main/notebooks/seismic.ipynb#Volume-Rendering)
* EarthCube2020: [citeable archive](https://www.essoar.org/doi/abs/10.1002/essoar.10506118.2) or the [repo](https://github.com/earthcube2020/ec20_havlin_etal).

both those rely on https://github.com/chrishavlin/yt_velmodel_vis

The EarthCube notebook ([direct link](https://nbviewer.jupyter.org/github/earthcube2020/ec20_havlin_etal/blob/master/notebook/ec20_havlin_etal.ipynb)) covers the interpolation above and goes into transfer functions in more detail.

### geo-volume rendering in practice:

New package!

https://github.com/chrishavlin/ytgeotools

streamlines all the above, with some additional analysis functionality.

import yt
from ytgeotools.seismology.datasets import XarrayGeoSpherical
filename = "IRIS/GYPSUM_percent.nc"

# to get a yt dataset for maps:
ds_yt = XarrayGeoSpherical(filename).load_uniform_grid()

# to get a yt dataset on interpolated cartesian grid:
ds = XarrayGeoSpherical(filename)
ds_yt_i = ds.interpolate_to_uniform_cartesian(
    ["dvs"],
    N=50,
    max_dist=50,
    return_yt=True,    
)

**should** be useable for any IRIS EMC netcdf file to get a yt dataset.

Installation not yet streamlined... need to install two packages from source in the following order:

1. https://github.com/yt-project/yt_idv
2. https://github.com/chrishavlin/ytgeotools

Need `cartopy`... easiest to use a conda environment.

Download/clone 1 and 2, make sure `conda` environment is active, then `cd` into `yt_idv` and install from source with

$ python pip install . 

Repeat for `ytgeotools`.

Test installation: from a python shell, try:

>>> import yt_idv
>>> import ytgeotools
>>> import yt

There may be some more missing depencies to install...

#### manual volume rendering

The AGU and EarthCube notebooks use the standard yt volume rendering interfaces. Excerpted from those notebooks, typical workflow to create an image would look like:

import yt
from ytgeotools.seismology.datasets import XarrayGeoSpherical

# load the cartesian dataset
ds_raw = XarrayGeoSpherical(filename)
ds = ds_raw.interpolate_to_uniform_cartesian(
    ["dvs"],
    N=50,
    max_dist=50,
    return_yt=True,    
)

# create the scene (loads full dataset into the scene for rendering)
sc = yt.create_scene(ds, "dvs") 

# adjust camera position and orientation (so "up" is surface)
x_c=np.mean(bbox[0])  
y_c=np.mean(bbox[1])
z_c=np.mean(bbox[2])
center_vec = np.array([x_c,y_c,z_c])
center_vec = center_vec / np.linalg.norm(center_vec)
pos=sc.camera.position
sc.camera.set_position(pos,north_vector=center_vec)

# adjust camera zoom 
zoom_factor=0.7 # < 1 zooms in
init_width=sc.camera.width
sc.camera.width = (init_width * zoom_factor)

# set transfer function
# initialize the tf object by setting the data bounds to consider
dvs_min=-8
dvs_max=8 
tf = yt.ColorTransferFunction((dvs_min,dvs_max))

# set gaussians to add 
TF_gaussians=[
    {'center':-.8,'width':.1,'RGBa':(1.,0.,0.,.5)},
    {'center':.5,'width':.2,'RGBa':(0.1,0.1,1.,.8)}
]

for gau in TF_gaussians:
    tf.add_gaussian(gau['center'],gau['width'],gau['RGBa'])
    
source = sc.sources['source_00']
source.set_transfer_function(tf)

# adjust resolution of rendering 
res = sc.camera.get_resolution()  
res_factor = 2
new_res = (int(res[0]*res_factor),int(res[1]*res_factor))
sc.camera.set_resolution(new_res)

# if in a notebook
sc.show()

# if in a script
sc.save()

A lot of trial and error for transfer functions, camera positioning. Make videos by adjusting the camera position, saving off the images and then stitching together with a different tool.

### interactive volume rendering

yt_idv!

import numpy as np
import yt_idv
from ytgeotools.seismology.datasets import XarrayGeoSpherical


def refill(vals):
    vals[np.isnan(vals)] = 0.0
    vals[vals > 0] = 0.0
    return vals


filename = "IRIS/NWUS11-S_percent.nc"
ds = XarrayGeoSpherical(filename)
ds_yt = ds.interpolate_to_uniform_cartesian(
    ["dvs"],
    N=50,
    max_dist=50,
    return_yt=True,
    rescale_coords=True, ## IMPORTANT!
    apply_functions=[refill, np.abs],
)

rc = yt_idv.render_context(height=800, width=800, gui=True)
sg = rc.add_scene(ds_yt, "dvs", no_ghost=True)
rc.run()

great for getting a sense of what's in the data!

caveats:
* max intensity and projection transfer functions only for now...
* need to install from source for now (maybe not later today?)

### annotations:

In manual plotting,

from yt.visualization.volume_rendering.api import LineSource, PointSource

Coordinates need to be in cartesian as well... The EarthCube/AGU notebooks have some helper functions to streamline this. Not yet in ytgeotools or yt_idv.

In manual plotting, adding the Point and Line sources can affect the overall brightness of the scene, often a lot of trial and error to adjust the opacity of the sources so that they do not wash out the volume rendering.