Properties of the cratons#
Requirement is to obtain some geophysical data about the cratonic lithosphere and compare to the same property for some other geological setting.
We can use the crust 1.0 geological regionalisation to query certain parts of the litho1.0 dataset and determine how lithospheric properties relate to the geological setting.
To do so, we need to be able to:
Construct a uniform discretisation of the globe for sampling
Obtain crustal type at those points
Select the points that meet certain criteria in the Crust 1.0 model
Extract properties from Litho1.0 for those points
Analyse
Plot
Stripy can produce the relevant sample points (we could, perhaps, use the original sample points for the litho1.0 model, or well- distributed points from a face-included icosahedron). Litho1.0 and Crust1.0 have the interpolation tools we need.
import litho1pt0 as litho1pt0
import stripy
import numpy as np
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 1
----> 1 import litho1pt0 as litho1pt0
2 import stripy
3 import numpy as np
ModuleNotFoundError: No module named 'litho1pt0'
even_mesh = stripy.spherical_meshes.icosahedral_mesh(include_face_points=True, refinement_levels=5)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 1
----> 1 even_mesh = stripy.spherical_meshes.icosahedral_mesh(include_face_points=True, refinement_levels=5)
NameError: name 'stripy' is not defined
number_of_mesh_points = even_mesh.npoints
print (number_of_mesh_points)
print (even_mesh.lats.shape)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 1
----> 1 number_of_mesh_points = even_mesh.npoints
2 print (number_of_mesh_points)
3 print (even_mesh.lats.shape)
NameError: name 'even_mesh' is not defined
latitudes_in_radians = even_mesh.lats
latitudes_in_degrees = np.degrees(latitudes_in_radians)
longitudes_in_radians = even_mesh.lons
longitudes_in_degrees = np.degrees(longitudes_in_radians)%360.0 - 180.0
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 1
----> 1 latitudes_in_radians = even_mesh.lats
2 latitudes_in_degrees = np.degrees(latitudes_in_radians)
4 longitudes_in_radians = even_mesh.lons
NameError: name 'even_mesh' is not defined
## Plot these points to check the mesh is reasonable
print (latitudes_in_degrees.min(), latitudes_in_degrees.max())
print (longitudes_in_degrees.min(), longitudes_in_degrees.max())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 3
1 ## Plot these points to check the mesh is reasonable
----> 3 print (latitudes_in_degrees.min(), latitudes_in_degrees.max())
4 print (longitudes_in_degrees.min(), longitudes_in_degrees.max())
NameError: name 'latitudes_in_degrees' is not defined
## Other version of this
# Make an empty array with the same size as the number of mesh points
crustype2 = np.empty(number_of_mesh_points)
crustype2 = np.empty_like(even_mesh.lats)
# Loop and fill the array with crust regionalisation values
for i in range(0,even_mesh.npoints):
if i%1000 == 0:
print("Iteration: {}".format(i))
crustype2[i] = litho1pt0.crust_type_at(lat=latitudes_in_degrees[i], lon=longitudes_in_degrees[i])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 4
1 ## Other version of this
2
3 # Make an empty array with the same size as the number of mesh points
----> 4 crustype2 = np.empty(number_of_mesh_points)
5 crustype2 = np.empty_like(even_mesh.lats)
7 # Loop and fill the array with crust regionalisation values
NameError: name 'np' is not defined
## Now need depth values from litho1pt0 geophysical model
litho1pt0.l1_layer_decode
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 3
1 ## Now need depth values from litho1pt0 geophysical model
----> 3 litho1pt0.l1_layer_decode
NameError: name 'litho1pt0' is not defined
l1 = litho1pt0.layer_depth(latitudes_in_degrees, longitudes_in_degrees, "LID-BOTTOM") *0.001
l2 = litho1pt0.layer_depth(latitudes_in_degrees, longitudes_in_degrees, "LID-TOP") *0.001
lthickness = (l1 - l2)
print("Litho thickness - {}, {}".format(lthickness.min(), lthickness.max()))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 1
----> 1 l1 = litho1pt0.layer_depth(latitudes_in_degrees, longitudes_in_degrees, "LID-BOTTOM") *0.001
2 l2 = litho1pt0.layer_depth(latitudes_in_degrees, longitudes_in_degrees, "LID-TOP") *0.001
4 lthickness = (l1 - l2)
NameError: name 'litho1pt0' is not defined
c1 = litho1pt0.layer_depth(latitudes_in_degrees, longitudes_in_degrees, "CRUST3-BOTTOM") *0.001
c2 = litho1pt0.layer_depth(latitudes_in_degrees, longitudes_in_degrees, "CRUST1-TOP") *0.001
s0 = litho1pt0.layer_depth(latitudes_in_degrees, longitudes_in_degrees, "SEDS1-TOP") *0.001
cthickness = (c1 - c2)
sthickness = (c2 - s0)
print("Crust thickness - {}, {}".format(cthickness.min(), cthickness.max()))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 1
----> 1 c1 = litho1pt0.layer_depth(latitudes_in_degrees, longitudes_in_degrees, "CRUST3-BOTTOM") *0.001
2 c2 = litho1pt0.layer_depth(latitudes_in_degrees, longitudes_in_degrees, "CRUST1-TOP") *0.001
4 s0 = litho1pt0.layer_depth(latitudes_in_degrees, longitudes_in_degrees, "SEDS1-TOP") *0.001
NameError: name 'litho1pt0' is not defined
## Is that sensible ?
%matplotlib inline
import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
global_extent = [-180.0, 180.0, -89, 89]
projection0 = ccrs.PlateCarree()
projection1 = ccrs.Orthographic(central_longitude=140.0, central_latitude=0.0, globe=None)
projection2 = ccrs.Mollweide()
projection3 = ccrs.Robinson()
base_projection = ccrs.PlateCarree()
fig = plt.figure(figsize=(12, 12), facecolor="none")
ax = plt.subplot(111, projection=projection0)
ax.set_global()
colormap = plt.cm.get_cmap('RdYlBu_r', 10)
m = ax.scatter(longitudes_in_degrees, latitudes_in_degrees, c=cthickness, cmap=colormap,
s=2.0, linewidth=0.0, transform=ccrs.PlateCarree())
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)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[10], line 3
1 ## Is that sensible ?
----> 3 get_ipython().run_line_magic('matplotlib', 'inline')
5 import cartopy
6 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'
crustype = litho1pt0.crust_type_at(lat=np.degrees(even_mesh.lats), lon=np.degrees(even_mesh.lons))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 1
----> 1 crustype = litho1pt0.crust_type_at(lat=np.degrees(even_mesh.lats), lon=np.degrees(even_mesh.lons))
NameError: name 'litho1pt0' is not defined
# litho1pt0.c1_region_descriptor in 2,3,4 are Archean terranes - this is an easy way to get those values
craton_cthickness = cthickness[np.logical_and(crustype<=4, crustype>=2)]
craton_lthickness = lthickness[np.logical_and(crustype<=4, crustype>=2)]
# litho1pt0.c1_region_descriptor in 5,6,7,8 are Proterozoic terranes - this is an easy way to get those values
proton_cthickness = cthickness[np.logical_and(crustype<=8, crustype>=5)]
proton_lthickness = lthickness[np.logical_and(crustype<=8, crustype>=5)]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 3
1 # litho1pt0.c1_region_descriptor in 2,3,4 are Archean terranes - this is an easy way to get those values
----> 3 craton_cthickness = cthickness[np.logical_and(crustype<=4, crustype>=2)]
4 craton_lthickness = lthickness[np.logical_and(crustype<=4, crustype>=2)]
6 # litho1pt0.c1_region_descriptor in 5,6,7,8 are Proterozoic terranes - this is an easy way to get those values
NameError: name 'cthickness' is not defined
print(craton_lthickness.mean())
print(proton_lthickness.mean())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 1
----> 1 print(craton_lthickness.mean())
2 print(proton_lthickness.mean())
NameError: name 'craton_lthickness' is not defined
for i, crustName in enumerate(litho1pt0.c1_region_descriptor):
if "Archean" in crustName:
print(i,"-",crustName)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 for i, crustName in enumerate(litho1pt0.c1_region_descriptor):
2 if "Archean" in crustName:
3 print(i,"-",crustName)
NameError: name 'litho1pt0' is not defined
Another worked example: Depth Age Relationship for the Oceans