Introduction to the stripy/python wrapper for litho 1.0#
Litho 1.0 is a global model of lithospheric properties that builds upon Crust 1.0. The original model was computed on an icosohedral triangulation (40962 points) and supplied with a C++ code to retrieve model information at any single lon, lat, depth location. To plot a map using Litho 1.0 data, numerous queries of the model are required for which repeated calls to the C++ code can be inefficient.
Using stripy to recreate the triangulation and provide interpolation through a (vectorised) numpy interface makes these queries straightforward. We provide these examples along with the
litho1pt0 source code to demonstrate a specific applicationo of the stripy wrapper.
Format for litho 1.0#
Litho 1.0 is supplied as a series of depths at each of the points of the triangulation. These layers represent distinct components of the lithospheric column that may not
be present everywhere. For efficiency of implementation, we have pre-processed the model to include all layers at every point but with zero thickness where any given layer is absent. This allows us to create global maps of layer depth or thickness using numpy arrays.
Below the lithosphere, the model blends seamlessly into a 1d Reference model (REF) which we do not include in this interface as it requires none of the triangulation and interpolation
capability of stripy.
References#
Pasyanos, M. E., T. G. Masters, G. Laske, and Z. Ma (2014), LITHO1.0: An updated crust and lithospheric model of the Earth, Journal of Geophysical Research-Solid Earth, 119(3), 2153–2173, doi:10.1002/2013JB010626.
Laske, G., G. Masters, and Z. Ma (2013), Update on CRUST1. 0—A 1-degree global model of Earth’s crust, Geophys Research Abstracts, 15, EGU2013–2658.
import litho1pt0 as litho
from pprint import pprint
import numpy as np
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 1
----> 1 import litho1pt0 as litho
2 from pprint import pprint
3 import numpy as np
ModuleNotFoundError: No module named 'litho1pt0'
litho._interpolator.npoints
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 1
----> 1 litho._interpolator.npoints
NameError: name 'litho' is not defined
print(" Layer keys \n")
pprint( litho.l1_layer_decode.items() )
print("\n Value keys \n")
pprint( litho.l1_data_decode.items() )
Layer keys
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 2
1 print(" Layer keys \n")
----> 2 pprint( litho.l1_layer_decode.items() )
3 print("\n Value keys \n")
4 pprint( litho.l1_data_decode.items() )
NameError: name 'pprint' is not defined
lats = np.array([0,10,20])
lons = np.array([0,0,0])
litho.layer_depth(lats,lons, layerID='ASTHENO-TOP')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 1
----> 1 lats = np.array([0,10,20])
2 lons = np.array([0,0,0])
4 litho.layer_depth(lats,lons, layerID='ASTHENO-TOP')
NameError: name 'np' is not defined
## make a global raster of some quantity
lonv, latv = np.meshgrid(np.linspace(-180,180,720), np.linspace(-89,89,356), sparse=False, indexing='xy')
l1 = litho.layer_depth(latv, lonv, "LID-BOTTOM")
l2 = litho.layer_depth(latv, lonv, "LID-TOP")
lthickness = (l1 - l2)*0.001
lab_depth = l1*0.001
l1 = litho.layer_depth(latv, lonv, "CRUST3-BOTTOM")
l2 = litho.layer_depth(latv, lonv, "CRUST1-TOP")
cthickness = (l1 - l2)*0.001
l1 = litho.layer_depth(latv, lonv, "LID-BOTTOM")
l2 = litho.layer_depth(latv, lonv, "CRUST1-TOP")
llthickness = (l1 - l2)*0.001
topo = litho.layer_depth(latv, lonv, "WATER-BOTTOM")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 3
1 ## make a global raster of some quantity
----> 3 lonv, latv = np.meshgrid(np.linspace(-180,180,720), np.linspace(-89,89,356), sparse=False, indexing='xy')
5 l1 = litho.layer_depth(latv, lonv, "LID-BOTTOM")
6 l2 = litho.layer_depth(latv, lonv, "LID-TOP")
NameError: name 'np' is not defined
%matplotlib inline
import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
global_extent = [-180.0, 180.0, -89, 89]
projection1 = ccrs.Orthographic(central_longitude=140.0, central_latitude=0.0, globe=None)
projection2 = ccrs.Mollweide()
projection3 = ccrs.Robinson()
base_projection = ccrs.PlateCarree()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[6], line 1
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')
3 import cartopy
4 import cartopy.crs as ccrs
File ~/.local/lib/python3.12/site-packages/IPython/core/interactiveshell.py:2482, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
2480 kwargs['local_ns'] = self.get_local_scope(stack_depth)
2481 with self.builtin_trap:
-> 2482 result = fn(*args, **kwargs)
2484 # The code below prevents the output from being displayed
2485 # when using magics with decorator @output_can_be_silenced
2486 # when the last Python token in the expression is a ';'.
2487 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):
File ~/.local/lib/python3.12/site-packages/IPython/core/magics/pylab.py:103, in PylabMagics.matplotlib(self, line)
98 print(
99 "Available matplotlib backends: %s"
100 % _list_matplotlib_backends_and_gui_loops()
101 )
102 else:
--> 103 gui, backend = self.shell.enable_matplotlib(args.gui)
104 self._show_matplotlib_backend(args.gui, backend)
File ~/.local/lib/python3.12/site-packages/IPython/core/interactiveshell.py:3667, in InteractiveShell.enable_matplotlib(self, gui)
3664 import matplotlib_inline.backend_inline
3666 from IPython.core import pylabtools as pt
-> 3667 gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)
3669 if gui != None:
3670 # If we have our first gui selection, store it
3671 if self.pylab_gui_select is None:
File ~/.local/lib/python3.12/site-packages/IPython/core/pylabtools.py:338, in find_gui_and_backend(gui, gui_select)
321 def find_gui_and_backend(gui=None, gui_select=None):
322 """Given a gui string return the gui and mpl backend.
323
324 Parameters
(...)
335 'WXAgg','Qt4Agg','module://matplotlib_inline.backend_inline','agg').
336 """
--> 338 import matplotlib
340 if _matplotlib_manages_backends():
341 backend_registry = matplotlib.backends.registry.backend_registry
ModuleNotFoundError: No module named 'matplotlib'
## Relative thickness of crust
fig = plt.figure(figsize=(12, 12), facecolor="none")
ax = plt.subplot(111, projection=projection2)
ax.set_global()
colormap = plt.cm.get_cmap('RdYlBu_r', 10)
"""
Possible values are:
Accent, Accent_r, Blues, Blues_r, BrBG, BrBG_r, BuGn, BuGn_r, BuPu, BuPu_r, CMRmap, CMRmap_r, Dark2, Dark2_r, GnBu, GnBu_r,
Greens, Greens_r, Greys, Greys_r, OrRd, OrRd_r, Oranges, Oranges_r, PRGn, PRGn_r, Paired, Paired_r, Pastel1, Pastel1_r, Pastel2,
Pastel2_r, PiYG, PiYG_r, PuBu, PuBuGn, PuBuGn_r, PuBu_r, PuOr, PuOr_r, PuRd, PuRd_r, Purples, Purples_r, RdBu, RdBu_r, RdGy, RdGy_r,
RdPu, RdPu_r, RdYlBu, RdYlBu_r, RdYlGn, RdYlGn_r, Reds, Reds_r, Set1, Set1_r, Set2, Set2_r, Set3, Set3_r, Spectral, Spectral_r, Vega10,
Vega10_r, Vega20, Vega20_r, Vega20b, Vega20b_r, Vega20c, Vega20c_r, Wistia, Wistia_r, YlGn, YlGnBu, YlGnBu_r, YlGn_r, YlOrBr, YlOrBr_r,
YlOrRd, YlOrRd_r, afmhot, afmhot_r, autumn, autumn_r, binary, binary_r, bone, bone_r, brg, brg_r, bwr, bwr_r, cool, cool_r, coolwarm,
coolwarm_r, copper, copper_r, cubehelix, cubehelix_r, flag, flag_r, gist_earth, gist_earth_r, gist_gray, gist_gray_r, gist_heat, gist_heat_r,
gist_ncar, gist_ncar_r, gist_rainbow, gist_rainbow_r, gist_stern, gist_stern_r, gist_yarg, gist_yarg_r, gnuplot, gnuplot2, gnuplot2_r,
gnuplot_r, gray, gray_r, hot, hot_r, hsv, hsv_r, inferno, inferno_r, jet, jet_r, magma, magma_r, nipy_spectral, nipy_spectral_r, ocean,
ocean_r, pink, pink_r, plasma, plasma_r, prism, prism_r, rainbow, rainbow_r, seismic, seismic_r, spectral, spectral_r, spring, spring_r,
summer, summer_r, tab10, tab10_r, tab20, tab20_r, tab20b, tab20b_r, tab20c, tab20c_r, terrain, terrain_r, viridis, viridis_r, winter, winter_r
"""
m = ax.imshow(cthickness/(cthickness+lthickness), origin='lower',
transform=base_projection,
extent=global_extent,
zorder=0,
cmap=colormap,
interpolation="gaussian")
plt.colorbar(mappable=m, orientation="horizontal", shrink=0.5)
ax.contour(lab_depth, origin='lower', levels=[250],
extent=global_extent, transform=base_projection, colors="#333333", alpha=0.75)
ax.add_feature(cartopy.feature.OCEAN, alpha=0.5, zorder=99, facecolor="#BBBBBB")
ax.coastlines(resolution="50m", zorder=100, linewidth=0.5)
# fig.savefig("RelativeCrustalThickness.png", dpi=600)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 3
1 ## Relative thickness of crust
----> 3 fig = plt.figure(figsize=(12, 12), facecolor="none")
4 ax = plt.subplot(111, projection=projection2)
5 ax.set_global()
NameError: name 'plt' is not defined
## 2nd Plot
fig = plt.figure(figsize=(12, 12), facecolor="none")
ax = plt.subplot(111, projection=projection2)
ax.set_global()
colormap = plt.cm.get_cmap('autumn', 10 )
m2 = ax.imshow(cthickness, origin='lower', transform=base_projection,
extent=global_extent, zorder=0, cmap=colormap, vmin=5, vmax=70,
interpolation="gaussian")
# m = ax.contourf(cthickness, origin='lower', levels=,
# cmap=colormap,
# extent=global_extent, transform=base_projection,
# extend="max")
cb = plt.colorbar(mappable=m2, orientation="horizontal", shrink=0.5)
# m = ax.contourf(lab_depth, origin='lower', levels=[250, 400],
# colors=["#FFFFFF"],
# extent=global_extent, transform=base_projection, extend="max",
# alpha=0.25)
m = ax.contour(lab_depth, origin='lower', levels=[250, ],
colors=["#FFFFFF"],
extent=global_extent, transform=base_projection, extend="max",
alpha=1.)
# ax.add_feature(cartopy.feature.OCEAN, alpha=0.5, zorder=99, facecolor="#BBBBBB")
ax.coastlines(resolution="50m", zorder=100, linewidth=0.5)
# plt.savefig("CrustalThickness.png", dpi=600 )
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 3
1 ## 2nd Plot
----> 3 fig = plt.figure(figsize=(12, 12), facecolor="none")
5 ax = plt.subplot(111, projection=projection2)
6 ax.set_global()
NameError: name 'plt' is not defined
## Contours of LAB
fig = plt.figure(figsize=(12, 12), facecolor="none")
ax = plt.subplot(111, projection=projection2)
ax.set_global()
m = ax.contourf(lab_depth, origin='lower', levels=[0, 130, 250, 400],
colors=[ "#FFFFFF", "#DDDDFF", "#9999FF"],
extent=global_extent, transform=base_projection, # extend="max",
linewidth=0.25)
# ax.contour(lab_depth, origin='lower', levels=[130, 250],
# extent=global_extent, transform=base_projection, colors="#555555")
plt.colorbar(mappable=m, orientation="horizontal", shrink=0.5)
ax.add_feature(cartopy.feature.OCEAN, alpha=0.5, zorder=99, facecolor="#BBBBBB")
ax.coastlines(resolution="50m", zorder=100, linewidth=0.5)
plt.savefig("LithosphereThickness.png", dpi=600)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 3
1 ## Contours of LAB
----> 3 fig = plt.figure(figsize=(12, 12), facecolor="none")
4 ax = plt.subplot(111, projection=projection2)
5 ax.set_global()
NameError: name 'plt' is not defined
Next example: Properties from Litho 1.0