This notebook demonstrates two ways of plotting seismic data from IRIS-EMC tomographic models:
The volume rendering here is based on the notebook presented at Havlin et al. 2020. In order to use the volume renderer, we'll need to interpolate the IRIS netcdf models from their geographic coordinates to cartesian coordinates, for which we'll use the https://github.com/chrishavlin/yt_velmodel_vis repository. The notebook from Havlin et al. 2020 describes the details of the interpolations and some additional annotations for volume rendering contained in the yt_velmodel_vis
repository, so for the present notebook, we've moved many of those function calls to the seismic_helper
module in the resources
directory for the present notebook. See the EarthCube notebook for a detailed explanation.
First, we load an IRIS EMC model into memory, use the yt_velmodel_vis.seis_model
class to interpolate to cartesian and then load a yt dataset using yt.load_uniform_grid
. We'll use shear wave velocity anomalies in the Western U.S. from James et al., 2011 (avaialable from IRIS here).
# imports and initialization
import os
import yt
import matplotlib.pyplot as plt
from yt_velmodel_vis import seis_model as SM, transferfunctions as TFs
from resources import seismic_helper as SH
# load interpolated data using the yt uniform grid loader (not editable)
# set the model file and the field to visualize
modelfile='NWUS11-S_percent.nc' # the model netCDF file
datafld='dvs' # the field to visualize, must be a variable in the netCDF file
# set the interpolation dictionary. If the interpolation for this model does
# not exist, SM.netcdf() will build it.
interp_dict={'field':datafld,'max_dist':50000,'res':[10000,10000,10000],
'input_units':'m','interpChunk':int(1e7)}
# load the model
model=SM.netcdf(modelfile,interp_dict)
# set some objects required for loading in yt
bbox = model.cart['bbox'] # the bounding box of interpolated cartesian grid
data={datafld:model.interp['data'][datafld]} # data container for yt scene
# load the data as a uniform grid, create the 3d scene
ds = yt.load_uniform_grid(data,data[datafld].shape,1.0,bbox=bbox,nprocs=1,
periodicity=(True,True,True),unit_system="mks")
yt : [INFO ] 2020-11-18 12:44:52,010 Parameters: current_time = 0.0 yt : [INFO ] 2020-11-18 12:44:52,011 Parameters: domain_dimensions = [287 271 237] yt : [INFO ] 2020-11-18 12:44:52,011 Parameters: domain_left_edge = [-5201282. -3445670.25 2740212.5 ] yt : [INFO ] 2020-11-18 12:44:52,012 Parameters: domain_right_edge = [-2328344.5 -735306.5625 5113796.5 ] yt : [INFO ] 2020-11-18 12:44:52,012 Parameters: cosmological_simulation = 0.0
Data loaded.
Volume rendering in yt requires that we specify a Scene
comprised of sources
representing the volumes of data that we wish to render, a transfer_function
that controls how a volume is sampled, and camera
settings that control orientation. The documentation on volume rendering here describes these components in more detail.
First, we set up a transfer function for the volume rendering following Havlin et al. 2020. This transfer function highlights a region of slow and fast anomalies, linearly varying the opacity and color within the slow and fast segments to emphasize interior structure:
# setting up transfer functions
tfOb = TFs.dv(data[datafld].ravel(),bounds=[-4,4])
# segment 1, slow anomalies
bnds=[-1.3,-.3]
TFseg=TFs.TFsegment(tfOb,bounds=bnds,cmap='OrRd_r')
alpha_o=0.95
Dalpha=-0.85
alpha=alpha_o + Dalpha/(bnds[1]-bnds[0]) * (TFseg.dvbins_c-bnds[0])
tfOb.addTFsegment(alpha,TFseg)
# segment 2, fast anomalies
bnds=[.1,.35]
TFseg=TFs.TFsegment(tfOb,bounds=bnds,cmap='winter_r')
alpha_o=.6
Dalpha=.4
alpha=alpha_o+ Dalpha/(bnds[1]-bnds[0]) * (TFseg.dvbins_c-bnds[0])
tfOb.addTFsegment(alpha,TFseg)
SH.plotTf(tfOb)
print("Ready to build scene")
Ready to build scene
Now we build our scene using the configure_scene
helper function, which sets the transfer function, camera orientation and adds annotations.
sc = SH.configure_scene(ds, datafld, model, bbox, tfOb.tf,res_factor=1.25)
/home/chavlin/miniconda3/envs/yt_dev/lib/python3.7/site-packages/geopandas/geodataframe.py:422: RuntimeWarning: Sequential read of iterator was interrupted. Resetting iterator. This can negatively impact the performance. for feature in features_lst:
Scene ready to render
The configure_scene
helper returns a scene that will render when we call sc.show()
:
sc.show(sigma_clip=1.5)
yt : [INFO ] 2020-11-18 12:45:32,123 Rendering scene (Can take a while). yt : [INFO ] 2020-11-18 12:45:32,268 Creating volume /home/chavlin/miniconda3/envs/yt_dev/lib/python3.7/site-packages/unyt/array.py:1649: RuntimeWarning: invalid value encountered in log10 out_arr = func(np.asarray(inp), out=out_func, **kwargs)
The annotations added by configure_scene
include state boundaries, plate boundaries and recent volcanos.
sc.save('./figures/seismic_render.png', sigma_clip = 1.5, render=False)
yt : [INFO ] 2020-11-18 15:14:11,754 Found previously rendered image to save. yt : [INFO ] 2020-11-18 15:14:11,757 Saving rendered image to ./figures/seismic_render.png
yt supports geographic transforms and plotting via Cartopy (link to yt docs). To take advantage of the geographic transforms, we can load the IRIS netcdf fiels directly into yt as uniform gridded data with the internal_geographic
geometry specification:
ds = yt.load_uniform_grid(data, sizes, 1.0, geometry=("internal_geographic", dims),bbox=bbox)
Before doing so, we first need to load our netcdf data into an in-memory flat data
dictionary, which we do here following the yt example notebook for using geographic projections.
import numpy as np
import yt
import os
import re
import cartopy.crs as ccrs
import cartopy.feature as cf
import netCDF4 as nc4
def get_data_path(arg):
if os.path.exists(arg):
return arg
else:
return os.path.join(yt.config.ytcfg.get("yt", "test_data_dir"), arg)
w_regex = re.compile(r'([a-zA-Z]+)(.*)')
def regex_parser(s):
try:
return "**".join(filter(None, w_regex.search(s).groups()))
except AttributeError:
return s
def get_dims(n):
dims = []
sizes = []
bbox = []
ndims = len(n.dimensions)
for dim in n.dimensions.keys():
size = n.variables[dim].size
if size > 1:
bbox.append([n.variables[dim][:].min(),
n.variables[dim][:].max()])
dims.append(n.variables[dim].name)
sizes.append(size)
bbox = np.array(bbox)
return dims, bbox, sizes, ndims
def prep_netcdf(n):
dims, bbox, sizes, ndims = get_dims(n)
data = {}
names = {}
for field, d in n.variables.items():
if d.ndim != ndims:
continue
units = n.variables[field].units
units = " * ".join(map(regex_parser, units.split()))
data[field] = (np.squeeze(d), str(units))
names[field] = n.variables[field].long_name.replace("_", " ")
return bbox, data, names, dims, sizes
At present, geographic plots in yt work best with global datasets, so for this example, we're using the global shear wave anomaly model of Simmons et al. 2010, available from IRIS here. So let's load and flatten our netcdf file:
fi = "GYPSUM_percent.nc"
n = nc4.Dataset(get_data_path(fi))
bbox, data, names, dims, sizes = prep_netcdf(n)
/home/chris/miniconda3/envs/yt_dev/lib/python3.7/site-packages/ipykernel_launcher.py:33: DeprecationWarning: tostring() is deprecated. Use tobytes() instead. /home/chris/miniconda3/envs/yt_dev/lib/python3.7/site-packages/ipykernel_launcher.py:36: DeprecationWarning: tostring() is deprecated. Use tobytes() instead.
Our global model covers the full range of latitude and longitude, but we have to currently set our bounds just under before creating our yt dataset:
bbox[1][0]=-89.9
bbox[1][1]=89.9
bbox[2][0]=-189.9
bbox[2][1]=189.9
ds = yt.load_uniform_grid(data, sizes, 1.0, geometry=("internal_geographic", dims),bbox=bbox)
yt : [INFO ] 2020-11-18 17:29:20,547 Parameters: current_time = 0.0 yt : [INFO ] 2020-11-18 17:29:20,547 Parameters: domain_dimensions = [100 181 361] yt : [INFO ] 2020-11-18 17:29:20,548 Parameters: domain_left_edge = [ 0. -89.90000153 -189.8999939 ] yt : [INFO ] 2020-11-18 17:29:20,548 Parameters: domain_right_edge = [ 2900. 89.90000153 189.8999939 ] yt : [INFO ] 2020-11-18 17:29:20,549 Parameters: cosmological_simulation = 0.0
Now that we have our dataset loaded, we can make a slice plot, setting our map projection using the set_mpl_projection
function:
target_depth = 100
center = ds.domain_center
p = yt.SlicePlot(ds,"depth",'dvs',center=[target_depth,center[1],center[2]]) # will be domain center
p.set_cmap('dvs', 'dusk_r')
p.set_zlim('dvs', -10, 10)
p.set_log('dvs', False)
p.annotate_title(f"Depth={target_depth} km")
p.set_mpl_projection(('PlateCarree', (), {'central_longitude':center[2]}))#{'central_longitude':-100, 'central_latitude':40}))
p._setup_plots()
p.plots['dvs'].axes.coastlines()
yt : [INFO ] 2020-11-18 17:22:29,623 Setting origin='native' for internal_geographic geometry. yt : [INFO ] 2020-11-18 17:22:29,630 xlim = -189.899994 189.899994 yt : [INFO ] 2020-11-18 17:22:29,631 ylim = -89.900002 89.900002 yt : [INFO ] 2020-11-18 17:22:29,634 xlim = -189.899994 189.899994 yt : [INFO ] 2020-11-18 17:22:29,635 ylim = -89.900002 89.900002 yt : [INFO ] 2020-11-18 17:22:29,638 Making a fixed resolution buffer of (('stream', 'dvs')) 800 by 800 yt : [WARNING ] 2020-11-18 17:22:29,731 Plot image for field ('stream', 'dvs') has both positive and negative values. Min = -7.125200, Max = 6.229640. yt : [WARNING ] 2020-11-18 17:22:29,732 Switching to symlog colorbar scaling unless linear scaling is specified later
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f7e8ad306d0>
p.show()
We can easily switch projections, in this case to a Mollweide projection:
p = yt.SlicePlot(ds,"depth",'dvs',center=[target_depth,center[1],center[2]]) # will be domain center
p.set_cmap('dvs', 'dusk_r')
p.set_zlim('dvs', -10, 10)
p.set_log('dvs', False)
p.annotate_title(f"Depth={target_depth} km")
p.set_mpl_projection(('Mollweide', (), {'central_longitude':center[2]}))#{'central_longitude':-100, 'central_latitude':40}))
p._setup_plots()
p.plots['dvs'].axes.coastlines()
yt : [INFO ] 2020-11-18 17:22:36,379 Setting origin='native' for internal_geographic geometry. yt : [INFO ] 2020-11-18 17:22:36,386 xlim = -189.899994 189.899994 yt : [INFO ] 2020-11-18 17:22:36,387 ylim = -89.900002 89.900002 yt : [INFO ] 2020-11-18 17:22:36,390 xlim = -189.899994 189.899994 yt : [INFO ] 2020-11-18 17:22:36,391 ylim = -89.900002 89.900002 yt : [INFO ] 2020-11-18 17:22:36,394 Making a fixed resolution buffer of (('stream', 'dvs')) 800 by 800 yt : [WARNING ] 2020-11-18 17:22:36,460 Plot image for field ('stream', 'dvs') has both positive and negative values. Min = -7.125200, Max = 6.229640. yt : [WARNING ] 2020-11-18 17:22:36,461 Switching to symlog colorbar scaling unless linear scaling is specified later
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f7e8a9ac210>
p.show()
yt allows access to the underlying Cartopy axes, so we can easily add data annotations directly to the axes using standard matplotlib/cartopy commands.
First, let's pull some data. Here, we're using obspy to query the IRIS catalogue for information about all earthquakes with a minimum magnitude of 5 from 2019-01-01 through 2020-11-01:
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
client = Client("IRIS")
starttime = UTCDateTime("2019-01-01")
endtime = UTCDateTime("2020-11-01")
cat = client.get_events(starttime=starttime, endtime=endtime,minmagnitude=5)
Now we can select the event locations:
lat_lon_dep = [[ev.origins[0].latitude,ev.origins[0].longitude,ev.origins[0].depth] for ev in cat.events]
lat_lon_dep = np.array(lat_lon_dep)
lat_lon_dep.shape
(2874, 3)
and now add them to our plot by first calling _setup_plots
and then accesing the plots
dictionary:
p = yt.SlicePlot(ds,"depth",'dvs',center=[target_depth,center[1],center[2]]) # will be domain center
p.set_cmap('dvs', 'dusk_r')
p.set_zlim('dvs', -10, 10)
p.set_log('dvs', False)
p.annotate_title(f"Depth={target_depth} km")
p.set_mpl_projection(('Mollweide', (), {'central_longitude':center[2]}))#{'central_longitude':-100, 'central_latitude':40}))
p._setup_plots()
p.plots['dvs'].axes.scatter(lat_lon_dep[:,1],lat_lon_dep[:,0],
transform=ccrs.PlateCarree(),color=[.8,.8,.8,.5])
p.plots['dvs'].axes.coastlines()
p.show()
yt : [INFO ] 2020-11-18 17:47:51,534 Setting origin='native' for internal_geographic geometry. yt : [INFO ] 2020-11-18 17:47:51,542 xlim = -189.899994 189.899994 yt : [INFO ] 2020-11-18 17:47:51,543 ylim = -89.900002 89.900002 yt : [INFO ] 2020-11-18 17:47:51,546 xlim = -189.899994 189.899994 yt : [INFO ] 2020-11-18 17:47:51,548 ylim = -89.900002 89.900002 yt : [INFO ] 2020-11-18 17:47:51,552 Making a fixed resolution buffer of (('stream', 'dvs')) 800 by 800 yt : [WARNING ] 2020-11-18 17:47:51,618 Plot image for field ('stream', 'dvs') has both positive and negative values. Min = -7.125200, Max = 6.229640. yt : [WARNING ] 2020-11-18 17:47:51,619 Switching to symlog colorbar scaling unless linear scaling is specified later
we can also use Cartopy's built in shapefile reader to add data.
In the following cell, we read the Harvard Global Volcanic Database shapefile (link) containing locations of Holocene volcanos:
from cartopy.io.shapereader import Reader
volc_file = get_data_path('harvard-glb-volc-shapefile/GLB_VOLC.shp')
volc_reader = Reader(volc_file)
which we can again add to our plot:
p = yt.SlicePlot(ds,"depth",'dvs',center=[target_depth,center[1],center[2]]) # will be domain center
p.set_cmap('dvs', 'dusk_r')
p.set_zlim('dvs', -10, 10)
p.set_log('dvs', False)
p.annotate_title(f"Depth={target_depth} km")
p.set_mpl_projection(('Mollweide', (), {'central_longitude':center[2]}))#{'central_longitude':-100, 'central_latitude':40}))
p._setup_plots()
points = list(volc_reader.geometries())
p.plots['dvs'].axes.scatter([point.x for point in points],
[point.y for point in points],
transform=ccrs.PlateCarree(),color='k')
p.plots['dvs'].axes.scatter(lat_lon_dep[:,1],lat_lon_dep[:,0],
transform=ccrs.PlateCarree(),color=[.8,.8,.8,.5])
p.plots['dvs'].axes.coastlines()
p.show()
yt : [INFO ] 2020-11-18 17:47:57,783 Setting origin='native' for internal_geographic geometry. yt : [INFO ] 2020-11-18 17:47:57,789 xlim = -189.899994 189.899994 yt : [INFO ] 2020-11-18 17:47:57,791 ylim = -89.900002 89.900002 yt : [INFO ] 2020-11-18 17:47:57,794 xlim = -189.899994 189.899994 yt : [INFO ] 2020-11-18 17:47:57,796 ylim = -89.900002 89.900002 yt : [INFO ] 2020-11-18 17:47:57,801 Making a fixed resolution buffer of (('stream', 'dvs')) 800 by 800 yt : [WARNING ] 2020-11-18 17:47:57,869 Plot image for field ('stream', 'dvs') has both positive and negative values. Min = -7.125200, Max = 6.229640. yt : [WARNING ] 2020-11-18 17:47:57,870 Switching to symlog colorbar scaling unless linear scaling is specified later
p.save('./figures/seismic_fixed_depth.png')
yt : [INFO ] 2020-11-18 17:48:01,008 Saving plot ./figures/seismic_fixed_depth.png
['./figures/seismic_fixed_depth.png']