--- 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: ```python= 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: ```python= cart_bbox = [ [min_x, max_x], [min_y, max_y], [min_z, max_z], ] ``` 3. create a uniform grid to cover bounding box: ```python= 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`: ```python 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. ```python 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: ```python= 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! ```python= 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, ```python= 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.
Last changed by  
 
chavlin
368

Read more

yt + Dask pYTEP

Table of Contents Overview Experiments in daskifying yt Development plan Some Final Notes yt and Dask: an overview In the past months, I've been investigating and working on integrating Dask into the yt codebase. This document provides an overview of my efforts to date but also is meant as a a preliminary YTEP (or pYTEP?) to solicit feedback from the yt community at an early stage before getting to far into the weeds of refactoring.

1/15/2021
yt + Dask particle IO

return to main post 1. (particle) data IO Full notebook available here. yt reads particle and grid-based data by iterating across the chunks, with frontend-specific IO functions. For gridded data, each frontend implements a _read_fluid_selection (e.g., yt.frontend.amrvac.AMRVACIOHandler._read_fluid_selection) that iterates over chunks and returns a flat dictionary with numpy arrays concatenated across each chunk. The particle data, frontends must implement a similar function, _read_particle_fields, that typically gets invoked within the BaseIOHanlder._read_particle_selection function. In both cases, the read functions accept the chunks iterator, the fields to read and a selector object: def _read_particle_fields(self, chunks, ptf, selector): def _read_fluid_selection(self, chunks, selector, fields, size):

1/15/2021
dask-unyt arrays

return to main post 3. dask-unyt arrays yt uses unyt to track and convert units so if we are using dask for IO and want to return delayed arrays, we need some level of dask-unyt support. In the notebook, working with unyt and dask, I demonstrate an initial prototype of a dask-unyt array. In this notebook, I create a custom dask collection by sublcassing the primary dask.array class and adding some unyt functionality in hidden sidecar attributes. This custom class is handled automatically by the dask scheduler, so that if we have a large dask array with a dask client running and we create our new dask-unyt array, e.g.: import dask.array as da from dask.distributed import Client client = Client(threads_per_worker=2, n_workers=2)

1/6/2021
yt + Dask profiles

return to main post 2. profile calculation: operations on chunked data The prototype dask-Gadget reader above returns the flat numpy arrays expected by _read_particle_selection(), but if instead we could return dask objects, we could easily build up parallel workflows. In the notebook here, I explored how to create a daskified-profile calculation under the assumption that our yt IO returns dask arrays. The notebook demonstrates two approaches. In the first, I use intrinsic dask array functions to recreate a yt profile using dask's implementation of np.digitize and then summing over the digitized bins. The performance isn't great, though, and it's not obvious how the approach could be extended to reproduce all of the functionality of ds.create_profile() in yt (multiple binning variables, returning multiple statistics at once, etc.). In the second approach, I instead work out how to apply the existing yt profile functions directly to the dask arrays. When you have a delayed dask array, you can loop over a chunk, apply a function to each chunk and then reduce the result across the chunks. In pseudo-code, creating a profile looks like:

1/6/2021
Read more from chavlin
Published on HackMD