yt_xarray and cartopy

This notebook demonstrates loading a geographic xarray dataset into yt. In addition to yt_xarray and xarray, this notebook requires cartopy and relies on the netcdf-gridded seismic tomography model of Schmand and Humphreys 2010, available from IRIS here.

import xarray as xr
import yt_xarray
import yt
from cartopy.feature import NaturalEarthFeature

yt_xarray provides a simple wrapper of the standard xarray open_dataset function that will check yt’s test_data_dir for data if the file is not found in the local path:

ds = yt_xarray.open_dataset("IRIS/wUS-SH-2010_percent.nc")

ds

<xarray.Dataset>
Dimensions:    (depth: 19, latitude: 93, longitude: 122)
Coordinates:
  * depth      (depth) float32 60.0 90.0 125.0 160.0 ... 700.0 760.0 820.0 885.0
  * latitude   (latitude) float32 27.5 27.75 28.0 28.25 ... 50.0 50.25 50.5
  * longitude  (longitude) float32 -125.8 -125.5 -125.2 ... -96.0 -95.75 -95.5
Data variables:
    dvp        (depth, latitude, longitude) float32 ...
    dvs        (depth, latitude, longitude) float32 ...
Attributes: (12/32)
    title:                         P and S teleseismic body-wave tomography o...
    id:                            wUS-SH-2010_percent
    summary:                       Teleseismic travel-time residuals from the...
    keywords:                      seismic, tomography, compressional wave, p...
    Conventions:                   CF-1.0
    Metadata_Conventions:          Unidata Dataset Discovery v1.0
    ...                            ...
    author_email:                  bschmandt@unm.edu
    author_institution:            Department of Earth and Planetary Science,...
    author_url:
    repository_name:               EMC
    repository_institution:        IRIS DMC
    repository_pid:                doi:10.17611/DP/9991760

xarray.Dataset

• Dimensions:
  • depth: 19
  • latitude: 93
  • longitude: 122
• Coordinates: (3)
  • depth
    (depth)
    float32
    60.0 90.0 125.0 ... 820.0 885.0

    long_name :

    depth below earth surface

    units :

    km

    positive :

    down

    array([ 60.,  90., 125., 160., 195., 230., 270., 310., 350., 395., 440., 490.,
           540., 590., 645., 700., 760., 820., 885.], dtype=float32)

  • latitude
    (latitude)
    float32
    27.5 27.75 28.0 ... 50.0 50.25 50.5

    long_name :

    Latitude; positive north

    units :

    degrees_north

    standard_name :

    latitude

    array([27.5 , 27.75, 28.  , 28.25, 28.5 , 28.75, 29.  , 29.25, 29.5 , 29.75,
           30.  , 30.25, 30.5 , 30.75, 31.  , 31.25, 31.5 , 31.75, 32.  , 32.25,
           32.5 , 32.75, 33.  , 33.25, 33.5 , 33.75, 34.  , 34.25, 34.5 , 34.75,
           35.  , 35.25, 35.5 , 35.75, 36.  , 36.25, 36.5 , 36.75, 37.  , 37.25,
           37.5 , 37.75, 38.  , 38.25, 38.5 , 38.75, 39.  , 39.25, 39.5 , 39.75,
           40.  , 40.25, 40.5 , 40.75, 41.  , 41.25, 41.5 , 41.75, 42.  , 42.25,
           42.5 , 42.75, 43.  , 43.25, 43.5 , 43.75, 44.  , 44.25, 44.5 , 44.75,
           45.  , 45.25, 45.5 , 45.75, 46.  , 46.25, 46.5 , 46.75, 47.  , 47.25,
           47.5 , 47.75, 48.  , 48.25, 48.5 , 48.75, 49.  , 49.25, 49.5 , 49.75,
           50.  , 50.25, 50.5 ], dtype=float32)

  • longitude
    (longitude)
    float32
    -125.8 -125.5 ... -95.75 -95.5

    long_name :

    Longitude; positive east

    units :

    degrees_east

    standard_name :

    longitude

    array([-125.75, -125.5 , -125.25, -125.  , -124.75, -124.5 , -124.25, -124.  ,
           -123.75, -123.5 , -123.25, -123.  , -122.75, -122.5 , -122.25, -122.  ,
           -121.75, -121.5 , -121.25, -121.  , -120.75, -120.5 , -120.25, -120.  ,
           -119.75, -119.5 , -119.25, -119.  , -118.75, -118.5 , -118.25, -118.  ,
           -117.75, -117.5 , -117.25, -117.  , -116.75, -116.5 , -116.25, -116.  ,
           -115.75, -115.5 , -115.25, -115.  , -114.75, -114.5 , -114.25, -114.  ,
           -113.75, -113.5 , -113.25, -113.  , -112.75, -112.5 , -112.25, -112.  ,
           -111.75, -111.5 , -111.25, -111.  , -110.75, -110.5 , -110.25, -110.  ,
           -109.75, -109.5 , -109.25, -109.  , -108.75, -108.5 , -108.25, -108.  ,
           -107.75, -107.5 , -107.25, -107.  , -106.75, -106.5 , -106.25, -106.  ,
           -105.75, -105.5 , -105.25, -105.  , -104.75, -104.5 , -104.25, -104.  ,
           -103.75, -103.5 , -103.25, -103.  , -102.75, -102.5 , -102.25, -102.  ,
           -101.75, -101.5 , -101.25, -101.  , -100.75, -100.5 , -100.25, -100.  ,
            -99.75,  -99.5 ,  -99.25,  -99.  ,  -98.75,  -98.5 ,  -98.25,  -98.  ,
            -97.75,  -97.5 ,  -97.25,  -97.  ,  -96.75,  -96.5 ,  -96.25,  -96.  ,
            -95.75,  -95.5 ], dtype=float32)

• Data variables: (2)
  • dvp
    (depth, latitude, longitude)
    float32
    ...

    long_name :

    Compressional Velocity perturbation

    display_name :

    dVp (%)

    units :

    percent

    [215574 values with dtype=float32]

  • dvs
    (depth, latitude, longitude)
    float32
    ...

    long_name :

    Shear Velocity perturbation

    display_name :

    dVs (%)

    units :

    percent

    [215574 values with dtype=float32]

• Attributes: (32)

  title :

  P and S teleseismic body-wave tomography of the western United States

  id :

  wUS-SH-2010_percent

  summary :

  Teleseismic travel-time residuals from the EarthScope Transportable Array and more than 1700 additional stations are inverted for 3-D velocity perturbations. The inversion uses frequency-dependent 3-D sensitivity kernels to map travel-time residuals, measured in multiple frequency bands, into velocity structure.

  keywords :

  seismic, tomography, compressional wave, p wave, shear wave, s wave, elsatic waveform

  Conventions :

  CF-1.0

  Metadata_Conventions :

  Unidata Dataset Discovery v1.0

  acknowledgment :

  Model was provided by , Brandon Schmandt, Department of Geological Sciences University of Oregon (now at Department of Earth and Planetary Sciences University of New Mexico)

  history :

  2019-12-05 IRIS DMC, updated geospatial fields, if needed, to float Wed Nov 20 10:29:11 2019: ncatted -O -a creator_email,global,m,c,product@iris.washington.edu wUS-SH-2010_percent.nc 2013-03-11 model converted to netCDF by IRIS EMC 2020-02-28 IRIS DMC, updated metadata to organize reference, author, repository and also add PID

  comment :

  model converted to netCDF by IRIS EMC

  geospatial_lat_min :

  27.5

  geospatial_lat_max :

  50.5

  geospatial_lat_units :

  degrees_north

  geospatial_lat_resolution :

  0.25

  geospatial_lon_min :

  -125.75

  geospatial_lon_max :

  -95.5

  geospatial_lon_units :

  degrees_east

  geospatial_lon_resolution :

  0.25

  geospatial_vertical_min :

  60

  geospatial_vertical_max :

  885

  geospatial_vertical_units :

  km

  geospatial_vertical_positive :

  down

  NCO :

  netCDF Operators version 4.7.5 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)

  netcdf_file :

  wUS-SH-2010_percent.nc

  reference :

  Schmandt and Humphreys (2010)

  reference_pid :

  doi:10.1016/j.epsl.2010.06.047

  author_name :

  Brandon Schmandt

  author_email :

  bschmandt@unm.edu

  author_institution :

  Department of Earth and Planetary Science, University of New Mexico

  author_url :

  repository_name :

  EMC

  repository_institution :

  IRIS DMC

  repository_pid :

  doi:10.17611/DP/9991760

to create a yt dataset that loads all the data variables, we simply call ds.yt.load_grid. For this dataset, our data is on a stretched grid, meaning the spacing in each dimension varies within a dimension (i.e., the depth spacing is changes), in which case we have to supply the use_callable=False flag to explicitly tell yt_xarray that we are OK with loading the data into memory (one of the current limitations of loading stretched grids):

yt_ds = ds.yt.load_grid(use_callable=False)

yt_xarray : [INFO ] 2023-02-06 15:56:02,810:  Inferred geometry type is geodetic. To override, use ds.yt.set_geometry
yt_xarray : [INFO ] 2023-02-06 15:56:02,816:  Attempting to detect if yt_xarray will require field interpolation:
yt_xarray : [INFO ] 2023-02-06 15:56:02,820:      stretched grid detected: yt_xarray will interpolate.
yt : [INFO     ] 2023-02-06 15:56:03,011 Parameters: current_time              = 0.0
yt : [INFO     ] 2023-02-06 15:56:03,012 Parameters: domain_dimensions         = [ 18  92 121]
yt : [INFO     ] 2023-02-06 15:56:03,014 Parameters: domain_left_edge          = [  60.     27.5  -125.75]
yt : [INFO     ] 2023-02-06 15:56:03,015 Parameters: domain_right_edge         = [885.   50.5 -95.5]
yt : [INFO     ] 2023-02-06 15:56:03,016 Parameters: cosmological_simulation   = 0

From here, yt will automatically use cartopy when plotting perpendicular to the vertical axis:

c = yt_ds.domain_center.copy()
c[0] = 150.
slc = yt.SlicePlot(yt_ds, "depth", ("stream", "dvs"), center = c)
slc.set_log(("stream", "dvs"), False)
slc.set_cmap(("stream", "dvs"), "magma_r")
slc.set_zlim(("stream", "dvs"), -8, 8)
slc._setup_plots()

states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
                             name='admin_1_states_provinces')
slc[("stream", "dvs")].axes.add_feature(states, edgecolor='gray')
slc.show()

yt : [INFO     ] 2023-02-06 15:56:03,137 xlim = -125.750000 -95.500000
yt : [INFO     ] 2023-02-06 15:56:03,137 ylim = 27.500000 50.500000
yt : [INFO     ] 2023-02-06 15:56:03,138 Setting origin='native' for internal_geographic geometry.
yt : [INFO     ] 2023-02-06 15:56:03,139 xlim = -125.750000 -95.500000
yt : [INFO     ] 2023-02-06 15:56:03,140 ylim = 27.500000 50.500000
yt : [INFO     ] 2023-02-06 15:56:03,149 Making a fixed resolution buffer of (('stream', 'dvs')) 800 by 800