Module gplately.oceans
A module to generate grids of seafloor age, seafloor spreading rate
and other oceanic data from the PlateReconstruction and
gplately.PlotToplogies objects.
Gridding methods in this module have been adapted from Simon Williams'
development repository for an
auto-age-gridding workflow, and are kept
within the SeafloorGrid object.
The sample jupyter notebook
10-SeafloorGrid
demonstrates how the functionalities within SeafloorGrid work. Below you can find
documentation for each of SeafloorGrid's functions.
SeafloorGrid Methodology
There are two main steps that SeafloorGrid follows to generate grids:
- Preparation for reconstruction by topologies
- Reconstruction by topologies
The preparation step involves building a:
- global domain of initial points that populate the seafloor at
max_time, - continental mask that separates ocean points from continent regions per timestep, and
- set of points that emerge to the left and right of mid-ocean ridge segments per timestep, as well as the z-value to allocate to these points.
First, the global domain of initial points is created using
stripy's
icosahedral triangulated mesh. The number of points in this mesh can be
controlled using a refinement_levels integer (the larger this integer,
the more resolved the continent masks will be).
These points are spatially partitioned by plate ID so they can be passed
into a
point-in-polygon routine.
This identifies points that lie within
continental polygon boundaries and those that are in the ocean. From this,
continental masks are built
per timestep, and the initial seed points are
allocated ages at the first reconstruction timestep max_time. Each point's
initial age is calculated by dividing its proximity to the nearest
MOR segment by half its assumed spreading rate. This spreading rate
(initial_ocean_mean_spreading_rate) is assumed to be uniform for all points.
These initial points momentarily fill the global ocean basin, and all have uniform spreading rates.
Thus, the spreading rate grid at max_time will be uniformly populated with the initial_ocean_mean_spreading_rate (mm/yr).
The age grid at max_time will look like a series of smooth, linear age gradients clearly partitioned by
tectonic plates with unique plate IDs:
Ridge "line" topologies
are resolved at each reconstruction time step and partitioned
into segments with a valid stage rotation. Each segment is further divided into points
at a specified ridge sampling spacing (ridge_sampling). Each point is
ascribed a latitude, longitude, spreading rate and age (from plate reconstruction
model files, as opposed to ages of the initial ocean mesh points), a point index
and the general z-value that will be gridded onto it.
Reconstruction by topologies involves determining which points are active and
inactive (collided with a continent or subducted at a trench) for each reconstruction
time step. This is done using a hidden object in PlateReconstruction called
ReconstructByTopologies.
If an ocean point with a certain velocity on one plate ID transitions into another
rigid plate ID at another timestep (with another velocity), the velocity difference
between both plates is calculated. The point may have subducted/collided with a continent
if this velocity difference is higher than a specified velocity threshold (which can be
controlled with subduction_collision_parameters). To ascertain whether the point
should be deactivated, a displacement test is conducted. If the proximity of the
point's previous time position to the polygon boundary it is approaching is higher than
a set distance threshold, then the point is far enough away from the boundary that it
cannot be subducted or consumed by it, and hence the point is still active. Otherwise,
it is deemed inactive and deleted from the ocean basin mesh.
With each reconstruction time step, points from mid-ocean ridges (which have more accurate spreading rates and attributed valid times) will spread across the ocean floor. Eventually, points will be pushed into continental boundaries or subduction zones, where they are deleted. Ideally, all initial ocean points (from the Stripy icosahedral mesh) should be deleted over time. However, not all will be deleted - such points typically reside near continental boundaries. This happens if the emerged ridge points do not spread far enough to "phase out" these points at collision regions - likely due to insufficient reconstruction detail. These undeleted points form artefacts of anomalously high seafloor age that append over the reconstruction time range.
Once reconstruction by topologies determines the ocean basin snapshot per timestep, a data frame of all longitudes, latitudes, seafloor ages, spreading rates and any other attributed z values will be written to a gridding input file per timestep.
Each active longitude, latitude and chosen z value (identified by a gridding input file
column index integer, i.e. 2 is seafloor age) is gridded using nearest-neighbour
interpolation and written to a netCDF4 format.
Classes
- SeafloorGrid
Expand source code
""" A module to generate grids of seafloor age, seafloor spreading rate
and other oceanic data from the `gplately.PlateReconstruction` and
`gplately.PlotToplogies` objects.
Gridding methods in this module have been adapted from Simon Williams'
development repository for an
[auto-age-gridding workflow](https://github.com/siwill22/agegrid-0.1), and are kept
within the `SeafloorGrid` object.
The sample jupyter notebook
[10-SeafloorGrid](https://github.com/GPlates/gplately/blob/master/Notebooks/10-SeafloorGrids.ipynb)
demonstrates how the functionalities within `SeafloorGrid` work. Below you can find
documentation for each of `SeafloorGrid`'s functions.
`SeafloorGrid` Methodology
--------------------------
There are two main steps that `SeafloorGrid` follows to generate grids:
1. Preparation for reconstruction by topologies
2. Reconstruction by topologies
The preparation step involves building a:
* global domain of initial points that populate the seafloor at `max_time`,
* continental mask that separates ocean points from continent regions per timestep, and
* set of points that emerge to the left and right of mid-ocean
ridge segments per timestep, as well as the z-value to allocate to these
points.
First, the global domain of initial points is created using
[stripy's](https://github.com/underworldcode/stripy/blob/master/stripy/spherical_meshes.py#L27)
icosahedral triangulated mesh. The number of points in this mesh can be
controlled using a `refinement_levels` integer (the larger this integer,
the more resolved the continent masks will be).
These points are spatially partitioned by plate ID so they can be passed
into a
[point-in-polygon routine](https://gplates.github.io/gplately/oceans.html#gplately.oceans.point_in_polygon_routine).
This identifies points that lie within
continental polygon boundaries and those that are in the ocean. From this,
[continental masks are built](https://gplates.github.io/gplately/oceans.html#gplately.oceans.SeafloorGrid.build_all_continental_masks)
per timestep, and the initial seed points are
allocated ages at the first reconstruction timestep `max_time`. Each point's
initial age is calculated by dividing its proximity to the nearest
MOR segment by half its assumed spreading rate. This spreading rate
(`initial_ocean_mean_spreading_rate`) is assumed to be uniform for all points.
These initial points momentarily fill the global ocean basin, and all have uniform spreading rates.
Thus, the spreading rate grid at `max_time` will be uniformly populated with the `initial_ocean_mean_spreading_rate` (mm/yr).
The age grid at `max_time` will look like a series of smooth, linear age gradients clearly partitioned by
tectonic plates with unique plate IDs:
[Ridge "line" topologies](https://gplates.github.io/gplately/oceans.html#gplately.oceans.SeafloorGrid.build_all_MOR_seedpoints)
are resolved at each reconstruction time step and partitioned
into segments with a valid stage rotation. Each segment is further divided into points
at a specified ridge sampling spacing (`ridge_sampling`). Each point is
ascribed a latitude, longitude, spreading rate and age (from plate reconstruction
model files, as opposed to ages of the initial ocean mesh points), a point index
and the general z-value that will be gridded onto it.
Reconstruction by topologies involves determining which points are active and
inactive (collided with a continent or subducted at a trench) for each reconstruction
time step. This is done using a hidden object in `PlateReconstruction` called
`ReconstructByTopologies`.
If an ocean point with a certain velocity on one plate ID transitions into another
rigid plate ID at another timestep (with another velocity), the velocity difference
between both plates is calculated. The point may have subducted/collided with a continent
if this velocity difference is higher than a specified velocity threshold (which can be
controlled with `subduction_collision_parameters`). To ascertain whether the point
should be deactivated, a displacement test is conducted. If the proximity of the
point's previous time position to the polygon boundary it is approaching is higher than
a set distance threshold, then the point is far enough away from the boundary that it
cannot be subducted or consumed by it, and hence the point is still active. Otherwise,
it is deemed inactive and deleted from the ocean basin mesh.
With each reconstruction time step, points from mid-ocean ridges (which have more
accurate spreading rates and attributed valid times) will spread across the ocean
floor. Eventually, points will be pushed into continental boundaries or subduction
zones, where they are deleted. Ideally, all initial ocean points (from the Stripy
icosahedral mesh) should be deleted over time. However, not all will be deleted -
such points typically reside near continental boundaries. This happens if the
emerged ridge points do not spread far enough to "phase out" these points at
collision regions - likely due to insufficient reconstruction detail. These
undeleted points form artefacts of anomalously high seafloor age that append
over the reconstruction time range.
Once reconstruction by topologies determines the ocean basin snapshot per timestep,
a data frame of all longitudes, latitudes, seafloor ages, spreading rates and any other
attributed z values will be written to a gridding input file per timestep.
Each active longitude, latitude and chosen z value (identified by a gridding input file
column index integer, i.e. `2` is seafloor age) is gridded using nearest-neighbour
interpolation and written to a netCDF4 format.
Classes
-------
* SeafloorGrid
"""
import glob
import logging
import os
import re
import warnings
import numpy as np
import pandas as pd
import pygplates
from . import grids, ptt, reconstruction, tools
from .ptt import separate_ridge_transform_segments
logger = logging.getLogger("gplately")
# -------------------------------------------------------------------------
# Auxiliary functions for SeafloorGrid
def create_icosahedral_mesh(refinement_levels):
"""Define a global point mesh with Stripy's
[icosahedral triangulated mesh](https://github.com/underworldcode/stripy/blob/294354c00dd72e085a018e69c345d9353c6fafef/stripy/spherical_meshes.py#L27)
and turn all mesh domains into pyGPlates MultiPointOnSphere types.
This global mesh will be masked with a set of continental or COB terrane
polygons to define the ocean basin at a given reconstruction time.
The `refinement_levels` integer is proportional to the resolution of the
mesh and the ocean/continent boundary.
Parameters
----------
refinement_levels : int
Refine the number of points in the triangulation. The larger the
refinement level, the sharper the ocean basin resolution.
Returns
-------
multi_point : instance of <pygplates.MultiPointOnSphere>
The longitues and latitudes that make up the icosahedral ocean mesh
collated into a MultiPointOnSphere object.
icosahedral_global_mesh : instance of <stripy.spherical_meshes.icosahedral_mesh>
The original global icosahedral triangulated mesh.
"""
import stripy
# Create the ocean basin mesh using Stripy's icosahedral spherical mesh
icosahedral_global_mesh = stripy.spherical_meshes.icosahedral_mesh(
refinement_levels, include_face_points=False, trisection=False, tree=False
)
# Get lons and lats of mesh, and turn them into a MultiPointOnSphere
lats_arr = np.rad2deg(icosahedral_global_mesh.lats)
lons_arr = np.rad2deg(icosahedral_global_mesh.lons)
multi_point = pygplates.MultiPointOnSphere(zip(lats_arr, lons_arr))
return multi_point, icosahedral_global_mesh
def ensure_polygon_geometry(reconstructed_polygons, rotation_model, time):
"""Ensure COB terrane/continental polygon geometries are polygons
with reconstruction plate IDs and valid times.
Notes
-----
This step must be done so that the initial set of ocean basin points
(the Stripy icosahedral mesh) can be partitioned into plates using
each reconstruction plate ID for the given plate `model`.
This allows for an oceanic point-in-continental
polygon query for every identified plate ID. See documentation for
`point_in_polygon_routine` for more details.
`ensure_polygon_geometry` works as follows:
COB terrane/continental polygons are assumed to have been reconstructed
already in `reconstructed_polygons` (a list of
type <pygplates.ReconstructedFeatureGeometry>). The list contents are
turned into a <pygplates.FeatureCollection> to be ascribed a
`PolygonOnSphere` geometry, a reconstruction plate ID, and a valid time.
Once finished, this feature collection is turned back into a list of
instance <pygplates.ReconstructedFeatureGeometry> and returned.
This revert must be completed for compatibility with the subsequent
point-in-polygon routine.
Parameters
----------
reconstructed_polygons : list of instance <pygplates.ReconstructedFeatureGeometry>
If used in `SeafloorGrid`, these are automatically obtained from the
`PlotTopologies.continents` attribute (the reconstructed continental
polygons at the current reconstruction time).
rotation_model : instance of <pygplates.RotationModel>
A parameter for turning the <pygplates.FeatureCollection> back into a
list of instance <pygplates.ReconstructedFeatureGeometry> for
compatibility with the point-in-polygon routine.
"""
continent_FeatCol = []
# self._PlotTopologies_object.continents
for n in reconstructed_polygons:
continent_FeatCol.append(n.get_feature())
polygon_feats = pygplates.FeatureCollection(continent_FeatCol)
# From GPRM's force_polygon_geometries(); set feature attributes
# like valid times and plate IDs to each masking polygon
polygons = []
for feature in polygon_feats:
for geom in feature.get_all_geometries():
polygon = pygplates.Feature(feature.get_feature_type())
polygon.set_geometry(pygplates.PolygonOnSphere(geom))
polygon.set_reconstruction_plate_id(feature.get_reconstruction_plate_id())
# Avoid features in COBTerranes with invalid time
if feature.get_valid_time()[0] >= feature.get_valid_time()[1]:
polygon.set_valid_time(
feature.get_valid_time()[0], feature.get_valid_time()[1]
)
polygons.append(polygon)
cobter_polygon_features = pygplates.FeatureCollection(polygons)
# Turn the feature collection back into ReconstructedFeatureGeometry
# objects otherwise it will not work with PIP
reconstructed_cobter_polygons = []
pygplates.reconstruct(
cobter_polygon_features, rotation_model, reconstructed_cobter_polygons, time
)
return reconstructed_cobter_polygons
def point_in_polygon_routine(multi_point, COB_polygons):
"""Perform Plate Tectonic Tools' point in polygon routine to partition
points in a `multi_point` MultiPointOnSphere feature based on whether
they are inside or outside the polygons in `COB_polygons`.
Notes
-----
Assuming the `COB_polygons` have passed through `ensure_polygon_geometry`,
each polygon should have a plate ID assigned to it.
This PIP routine serves two purposes for `SeafloorGrid`:
1) It identifies continental regions in the icosahedral global mesh
MultiPointOnSphere feature and 'erases' in-continent oceanic points
for the construction of a continental mask at each timestep;
2) It identifies oceanic points in the icosahedral global mesh.
These points will be passed to a function that calculates each point's
proximity to its nearest MOR segment (if any) within the polygonal domain
of its allocated plate ID. Each distance is divided by half the
`initial_ocean_mean_spreading_rate` (an attribute of `SeafloorGrids`) to
determine a simplified seafloor age for each point.
Number 2) only happens once at the start of the gridding process to
momentarily fill the gridding region with initial ocean points that have
set ages (albeit not from a plate model file). After multiple time steps
of reconstruction, the ocean basin will be filled with new points (with
plate-model prescribed ages) that emerge from ridge topologies.
Returns
-------
pygplates.MultiPointOnSphere(points_in_arr) : instance <pygplates.MultiPointOnSphere>
Point features that are within COB terrane polygons.
pygplates.MultiPointOnSphere(points_out_arr) : instance <pygplates.MultiPointOnSphere>
Point features that are outside COB terrane polygons.
zvals : list
A binary list. If an entry is == 0, its corresponing point in the
MultiPointOnSphere object is on the ocean. If == 1, the point is
in the COB terrane polygon.
"""
# Convert MultiPointOnSphere to array of PointOnSphere
multi_point = np.array(multi_point.get_points(), dtype="object")
# Collect reconstructed geometries of continental polygons
polygons = np.empty(len(COB_polygons), dtype="object")
for ind, i in enumerate(COB_polygons):
if isinstance(i, pygplates.ReconstructedFeatureGeometry):
geom = i.get_reconstructed_geometry()
elif isinstance(i, pygplates.GeometryOnSphere):
geom = i
else: # e.g. ndarray of coordinates
geom = pygplates.PolygonOnSphere(i)
polygons[ind] = geom
proxies = np.ones(polygons.size)
pip_result = ptt.utils.points_in_polygons.find_polygons(
multi_point, polygons, proxies, all_polygons=False
) # 1 for points in polygons, None for points outside
zvals = np.array(
pip_result,
dtype="float",
).ravel()
zvals[np.isnan(zvals)] = 0.0
zvals = zvals.astype("int")
points_in_arr = multi_point[zvals == 1]
points_out_arr = multi_point[zvals != 1]
return (
pygplates.MultiPointOnSphere(points_in_arr),
pygplates.MultiPointOnSphere(points_out_arr),
zvals,
)
def _deg2pixels(deg_res, deg_min, deg_max):
return int(np.floor((deg_max - deg_min) / deg_res)) + 1
def _pixels2deg(spacing_pixel, deg_min, deg_max):
return (deg_max - deg_min) / np.floor(int(spacing_pixel - 1))
class SeafloorGrid(object):
"""A class to generate grids that track data atop global ocean basin points
(which emerge from mid ocean ridges) through geological time.
Parameters
----------
PlateReconstruction_object : instance of <gplately.PlateReconstruction>
A GPlately PlateReconstruction object with a <pygplates.RotationModel> and
a <pygplates.FeatureCollection> containing topology features.
PlotTopologies_object : instance of <gplately.PlotTopologies>
A GPlately PlotTopologies object with a continental polygon or COB terrane
polygon file to mask grids with.
max_time : float
The maximum time for age gridding.
min_time : float
The minimum time for age gridding.
ridge_time_step : float
The delta time for resolving ridges (and thus age gridding).
save_directory : str, default None'
The top-level directory to save all outputs to.
file_collection : str, default None
A string to identify the plate model used (will be automated later).
refinement_levels : int, default 5
Control the number of points in the icosahedral mesh (higher integer
means higher resolution of continent masks).
ridge_sampling : float, default 0.5
Spatial resolution (in degrees) at which points that emerge from ridges are tessellated.
extent : list of float or int, default [-180.,180.,-90.,90.]
A list containing the mininum longitude, maximum longitude, minimum latitude and
maximum latitude extents for all masking and final grids.
grid_spacing : float, default None
The degree spacing/interval with which to space grid points across all masking and
final grids. If `grid_spacing` is provided, all grids will use it. If not,
`grid_spacing` defaults to 0.1.
subduction_collision_parameters : len-2 tuple of float, default (5.0, 10.0)
A 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary
per My in kms/my)
initial_ocean_mean_spreading_rate : float, default 75.
A spreading rate to uniformly allocate to points that define the initial ocean
basin. These points will have inaccurate ages, but most of them will be phased
out after points with plate-model prescribed ages emerge from ridges and spread
to push them towards collision boundaries (where they are deleted).
resume_from_checkpoints : bool, default False
If set to `True`, and the gridding preparation stage (continental masking and/or
ridge seed building) is interrupted, SeafloorGrids will resume gridding preparation
from the last successful preparation time.
If set to `False`, SeafloorGrids will automatically overwrite all files in
`save_directory` if re-run after interruption, or normally re-run, thus beginning
gridding preparation from scratch. `False` will be useful if data allocated to the
MOR seed points need to be augmented.
zval_names : list of str
A list containing string labels for the z values to attribute to points.
Will be used as column headers for z value point dataframes.
continent_mask_filename : str
An optional parameter pointing to the full path to a continental mask for each timestep.
Assuming the time is in the filename, i.e. "/path/to/continent_mask_0Ma.nc", it should be
passed as "/path/to/continent_mask_{}Ma.nc" with curly brackets. Include decimal formatting
if needed.
"""
def __init__(
self,
PlateReconstruction_object,
PlotTopologies_object,
max_time,
min_time,
ridge_time_step,
save_directory="agegrids",
file_collection=None,
refinement_levels=5,
ridge_sampling=0.5,
extent=(-180, 180, -90, 90),
grid_spacing=None,
subduction_collision_parameters=(5.0, 10.0),
initial_ocean_mean_spreading_rate=75.0,
resume_from_checkpoints=False,
zval_names=("SPREADING_RATE",),
continent_mask_filename=None,
):
# Provides a rotation model, topology features and reconstruction time for
# the SeafloorGrid
self.PlateReconstruction_object = PlateReconstruction_object
self.rotation_model = self.PlateReconstruction_object.rotation_model
self.topology_features = self.PlateReconstruction_object.topology_features
self._PlotTopologies_object = PlotTopologies_object
save_directory = str(save_directory)
if not os.path.isdir(save_directory):
logger.info(
"Output directory does not exist; creating now: " + save_directory
)
os.makedirs(save_directory, exist_ok=True)
self.save_directory = save_directory
if file_collection is not None:
file_collection = str(file_collection)
self.file_collection = file_collection
# Topological parameters
self.refinement_levels = refinement_levels
self.ridge_sampling = ridge_sampling
self.subduction_collision_parameters = subduction_collision_parameters
self.initial_ocean_mean_spreading_rate = initial_ocean_mean_spreading_rate
# Gridding parameters
self.extent = extent
# A list of degree spacings that allow an even division of the global lat-lon extent.
divisible_degree_spacings = [0.1, 0.25, 0.5, 0.75, 1.0]
if grid_spacing:
# If the provided degree spacing is in the list of permissible spacings, use it
# and prepare the number of pixels in x and y (spacingX and spacingY)
if grid_spacing in divisible_degree_spacings:
self.grid_spacing = grid_spacing
self.spacingX = _deg2pixels(
grid_spacing, self.extent[0], self.extent[1]
)
self.spacingY = _deg2pixels(
grid_spacing, self.extent[2], self.extent[3]
)
# If the provided spacing is >>1 degree, use 1 degree
elif grid_spacing >= divisible_degree_spacings[-1]:
self.grid_spacing = divisible_degree_spacings[-1]
self.spacingX = _deg2pixels(
divisible_degree_spacings[-1], self.extent[0], self.extent[1]
)
self.spacingY = _deg2pixels(
divisible_degree_spacings[-1], self.extent[2], self.extent[3]
)
with warnings.catch_warnings():
warnings.simplefilter("always")
warnings.warn(
f"The provided grid_spacing of {grid_spacing} is quite large. To preserve the grid resolution, a {self.grid_spacing} degree spacing has been employed instead."
)
# If the provided degree spacing is not in the list of permissible spacings, but below
# a degree, find the closest permissible degree spacing. Use this and find
# spacingX and spacingY.
else:
for divisible_degree_spacing in divisible_degree_spacings:
# The tolerance is half the difference between consecutive divisible spacings.
# Max is 1 degree for now - other integers work but may provide too little of a
# grid resolution.
if abs(grid_spacing - divisible_degree_spacing) <= 0.125:
new_deg_res = divisible_degree_spacing
self.grid_spacing = new_deg_res
self.spacingX = _deg2pixels(
new_deg_res, self.extent[0], self.extent[1]
)
self.spacingY = _deg2pixels(
new_deg_res, self.extent[2], self.extent[3]
)
with warnings.catch_warnings():
warnings.simplefilter("always")
warnings.warn(
f"The provided grid_spacing of {grid_spacing} does not cleanly divide into the global extent. A degree spacing of {self.grid_spacing} has been employed instead."
)
else:
# If a spacing degree is not provided, use default
# resolution and get default spacingX and spacingY
self.grid_spacing = 0.1
self.spacingX = 3601
self.spacingY = 1801
self.resume_from_checkpoints = resume_from_checkpoints
# Temporal parameters
self._max_time = float(max_time)
self.min_time = float(min_time)
self.ridge_time_step = float(ridge_time_step)
self.time_array = np.arange(
self._max_time, self.min_time - 0.1, -self.ridge_time_step
)
# If PlotTopologies' time attribute is not equal to the maximum time in the
# seafloor grid reconstruction tree, make it equal. This will ensure the time
# for continental masking is consistent.
if self._PlotTopologies_object.time != self._max_time:
self._PlotTopologies_object.time = self._max_time
# Essential features and meshes for the SeafloorGrid
self.continental_polygons = ensure_polygon_geometry(
self._PlotTopologies_object.continents, self.rotation_model, self._max_time
)
self._PlotTopologies_object.continents = PlotTopologies_object.continents
(
self.icosahedral_multi_point,
self.icosahedral_global_mesh,
) = create_icosahedral_mesh(self.refinement_levels)
# Z value parameters
self.zval_names = zval_names
self.default_column_headers = [
"CURRENT_LONGITUDES",
"CURRENT_LATITUDES",
"SEAFLOOR_AGE",
"BIRTH_LAT_SNAPSHOT",
"POINT_ID_SNAPSHOT",
]
self.total_column_headers = np.concatenate(
[self.default_column_headers, self.zval_names]
)
# Filename for continental masks that the user can provide instead of building it here
self.continent_mask_filename = continent_mask_filename
# If the user provides a continental mask filename, we need to downsize the mask
# resolution for when we create the initial ocean mesh. The mesh does not need to be high-res.
if self.continent_mask_filename is not None:
# Determine which percentage to use to scale the continent mask resolution at max time
def _map_res_to_node_percentage(self, continent_mask_filename):
maskY, maskX = grids.read_netcdf_grid(
continent_mask_filename.format(self._max_time)
).shape
mask_deg = _pixels2deg(maskX, self.extent[0], self.extent[1])
if mask_deg <= 0.1:
percentage = 0.1
elif mask_deg <= 0.25:
percentage = 0.3
elif mask_deg <= 0.5:
percentage = 0.5
elif mask_deg < 0.75:
percentage = 0.6
elif mask_deg >= 1:
percentage = 0.75
return mask_deg, percentage
_, self.percentage = _map_res_to_node_percentage(
self, self.continent_mask_filename
)
# Allow SeafloorGrid time to be updated, and to update the internally-used
# PlotTopologies' time attribute too. If PlotTopologies is used outside the
# object, its `time` attribute is not updated.
@property
def max_time(self):
"""The reconstruction time."""
return self._max_time
@property
def PlotTopologiesTime(self):
return self._PlotTopologies_object.time
@max_time.setter
def max_time(self, var):
if var >= 0:
self.update_time(var)
else:
raise ValueError("Enter a valid time >= 0")
def update_time(self, max_time):
self._max_time = float(max_time)
self._PlotTopologies_object.time = float(max_time)
def _collect_point_data_in_dataframe(
self, pygplates_featurecollection, zval_ndarray, time
):
"""At a given timestep, create a pandas dataframe holding all attributes of point features.
Rather than store z values as shapefile attributes, store them in a dataframe indexed by
feature ID.
"""
# Turn the zval_ndarray into a numPy array
zval_ndarray = np.array(zval_ndarray)
feature_id = []
for feature in pygplates_featurecollection:
feature_id.append(str(feature.get_feature_id()))
# Prepare the zval ndarray (can be of any shape) to be saved with default point data
zvals_to_store = {}
# If only one zvalue (fow now, spreading rate)
if zval_ndarray.ndim == 1:
zvals_to_store[self.zval_names[0]] = zval_ndarray
data_to_store = [zvals_to_store[i] for i in zvals_to_store]
else:
for i in zval_ndarray.shape[1]:
zvals_to_store[self.zval_names[i]] = [
list(j) for j in zip(*zval_ndarray)
][i]
data_to_store = [zvals_to_store[i] for i in zvals_to_store]
basename = "point_data_dataframe_{}Ma".format(time)
if self.file_collection is not None:
basename = "{}_{}".format(self.file_collection, basename)
filename = os.path.join(self.save_directory, basename)
np.savez_compressed(filename, FEATURE_ID=feature_id, *data_to_store)
return
def create_initial_ocean_seed_points(self):
"""Create the initial ocean basin seed point domain (at `max_time` only)
using Stripy's icosahedral triangulation with the specified
`self.refinement_levels`.
The ocean mesh starts off as a global-spanning Stripy icosahedral mesh.
`create_initial_ocean_seed_points` passes the automatically-resolved-to-current-time
continental polygons from the `PlotTopologies_object`'s `continents` attribute
(which can be from a COB terrane file or a continental polygon file) into
Plate Tectonic Tools' point-in-polygon routine. It identifies ocean basin points
that lie:
* outside the polygons (for the ocean basin point domain)
* inside the polygons (for the continental mask)
Points from the mesh outside the continental polygons make up the ocean basin seed
point mesh. The masked mesh is outputted as a compressed GPML (GPMLZ) file with
the filename: "ocean_basin_seed_points_{}Ma.gpmlz" if a `save_directory` is passed.
Otherwise, the mesh is returned as a pyGPlates FeatureCollection object.
Notes
-----
This point mesh represents ocean basin seafloor that was produced
before `SeafloorGrid.max_time`, and thus has unknown properties like valid
time and spreading rate. As time passes, the plate reconstruction model sees
points emerging from MORs. These new points spread to occupy the ocean basins,
moving the initial filler points closer to subduction zones and continental
polygons with which they can collide. If a collision is detected by
`PlateReconstruction`s `ReconstructByTopologies` object, these points are deleted.
Ideally, if a reconstruction tree spans a large time range, **all** initial mesh
points would collide with a continent or be subducted, leaving behind a mesh of
well-defined MOR-emerged ocean basin points that data can be attributed to.
However, some of these initial points situated close to contiental boundaries are
retained through time - these form point artefacts with anomalously high ages. Even
deep-time plate models (e.g. 1 Ga) will have these artefacts - removing them would
require more detail to be added to the reconstruction model.
Returns
-------
ocean_basin_point_mesh : pygplates.FeatureCollection of pygplates.MultiPointOnSphere
A feature collection of point objects on the ocean basin.
"""
if self.continent_mask_filename is None:
# Ensure COB terranes at max time have reconstruction IDs and valid times
COB_polygons = ensure_polygon_geometry(
self._PlotTopologies_object.continents,
self.rotation_model,
self._max_time,
)
# zval is a binary array encoding whether a point
# coordinate is within a COB terrane polygon or not.
# Use the icosahedral mesh MultiPointOnSphere attribute
_, ocean_basin_point_mesh, zvals = point_in_polygon_routine(
self.icosahedral_multi_point, COB_polygons
)
# Plates to partition with
plate_partitioner = pygplates.PlatePartitioner(
COB_polygons,
self.rotation_model,
)
# Plate partition the ocean basin points
meshnode_feature = pygplates.Feature(
pygplates.FeatureType.create_from_qualified_string("gpml:MeshNode")
)
meshnode_feature.set_geometry(
ocean_basin_point_mesh
# multi_point
)
ocean_basin_meshnode = pygplates.FeatureCollection(meshnode_feature)
paleogeography = plate_partitioner.partition_features(
ocean_basin_meshnode,
partition_return=pygplates.PartitionReturn.separate_partitioned_and_unpartitioned,
properties_to_copy=[pygplates.PropertyName.gpml_shapefile_attributes],
)
ocean_points = paleogeography[1] # Separate those inside polygons
continent_points = paleogeography[0] # Separate those outside polygons
# If a set of continent masks was passed, we can use max_time's continental
# mask to build the initial profile of seafloor age.
else:
max_time_cont_mask = grids.Raster(
self.continent_mask_filename.format(self._max_time)
)
# If the input grid is at 0.5 degree uniform spacing, then the input
# grid is 7x more populated than a 6-level stripy icosahedral mesh and
# using this resolution for the initial ocean mesh will dramatically slow down
# reconstruction by topologies.
# Scale down the resolution based on the input mask resolution
# (percentage was found in __init__.)
max_time_cont_mask.resize(
int(max_time_cont_mask.shape[0] * self.percentage),
int(max_time_cont_mask.shape[1] * self.percentage),
inplace=True,
)
lat = np.linspace(-90, 90, max_time_cont_mask.shape[0])
lon = np.linspace(-180, 180, max_time_cont_mask.shape[1])
llon, llat = np.meshgrid(lon, lat)
mask_inds = np.where(max_time_cont_mask.data.flatten() == 0)
mask_vals = max_time_cont_mask.data.flatten()
mask_lon = llon.flatten()[mask_inds]
mask_lat = llat.flatten()[mask_inds]
ocean_pt_feature = pygplates.Feature()
ocean_pt_feature.set_geometry(
pygplates.MultiPointOnSphere(zip(mask_lat, mask_lon))
)
ocean_points = [ocean_pt_feature]
# Now that we have ocean points...
# Determine age of ocean basin points using their proximity to MOR features
# and an assumed globally-uniform ocean basin mean spreading rate.
# We need resolved topologies at the `max_time` to pass into the proximity
# function
resolved_topologies = []
shared_boundary_sections = []
pygplates.resolve_topologies(
self.topology_features,
self.rotation_model,
resolved_topologies,
self._max_time,
shared_boundary_sections,
)
pX, pY, pZ = tools.find_distance_to_nearest_ridge(
resolved_topologies,
shared_boundary_sections,
ocean_points,
)
# Divide spreading rate by 2 to use half the mean spreading rate
pAge = np.array(pZ) / (self.initial_ocean_mean_spreading_rate / 2.0)
initial_ocean_point_features = []
initial_ocean_multipoints = []
for point in zip(pX, pY, pAge):
point_feature = pygplates.Feature()
point_feature.set_geometry(pygplates.PointOnSphere(point[1], point[0]))
# Add 'time' to the age at the time of computation, to get the valid time in Ma
point_feature.set_valid_time(point[2] + self._max_time, -1)
# For now: custom zvals are added as shapefile attributes - will attempt pandas data frames
# point_feature = set_shapefile_attribute(point_feature, self.initial_ocean_mean_spreading_rate, "SPREADING_RATE") # Seems like static data
initial_ocean_point_features.append(point_feature)
initial_ocean_multipoints.append(point_feature.get_geometry())
# print(initial_ocean_point_features)
multi_point_feature = pygplates.MultiPointOnSphere(initial_ocean_multipoints)
basename = "ocean_basin_seed_points_{}_RLs_{}Ma.gpmlz".format(
self.refinement_levels,
self._max_time,
)
if self.file_collection is not None:
basename = "{}_{}".format(self.file_collection, basename)
output_filename = os.path.join(self.save_directory, basename)
initial_ocean_feature_collection = pygplates.FeatureCollection(
initial_ocean_point_features
)
initial_ocean_feature_collection.write(output_filename)
# Collect all point feature data into a pandas dataframe
self._collect_point_data_in_dataframe(
initial_ocean_feature_collection,
np.array(
[self.initial_ocean_mean_spreading_rate] * len(pX)
), # for now, spreading rate is one zvalue for initial ocean points. will other zvalues need to have a generalised workflow?
self._max_time,
)
return (
pygplates.FeatureCollection(initial_ocean_point_features),
multi_point_feature,
)
def _get_mid_ocean_ridge_seedpoints(self, time_array):
# Topology features from `PlotTopologies`.
topology_features_extracted = pygplates.FeaturesFunctionArgument(
self.topology_features
)
# Create a mask for each timestep
if time_array[0] != self._max_time:
print(
"MOR seed point building interrupted - resuming at {} Ma!".format(
time_array[0]
)
)
for time in time_array:
# Points and their z values that emerge from MORs at this time.
shifted_mor_points = []
point_spreading_rates = []
# Resolve topologies to the current time.
resolved_topologies = []
shared_boundary_sections = []
pygplates.resolve_topologies(
topology_features_extracted.get_features(),
self.rotation_model,
resolved_topologies,
time,
shared_boundary_sections,
)
# pygplates.ResolvedTopologicalSection objects.
for shared_boundary_section in shared_boundary_sections:
if (
shared_boundary_section.get_feature().get_feature_type()
== pygplates.FeatureType.create_gpml("MidOceanRidge")
):
spreading_feature = shared_boundary_section.get_feature()
# Find the stage rotation of the spreading feature in the
# frame of reference of its geometry at the current
# reconstruction time (the MOR is currently actively spreading).
# The stage pole can then be directly geometrically compared
# to the *reconstructed* spreading geometry.
stage_rotation = separate_ridge_transform_segments.get_stage_rotation_for_reconstructed_geometry(
spreading_feature, self.rotation_model, time
)
if not stage_rotation:
# Skip current feature - it's not a spreading feature.
continue
# Get the stage pole of the stage rotation.
# Note that the stage rotation is already in frame of
# reference of the *reconstructed* geometry at the spreading time.
stage_pole, _ = stage_rotation.get_euler_pole_and_angle()
# One way rotates left and the other right, but don't know
# which - doesn't matter in our example though.
rotate_slightly_off_mor_one_way = pygplates.FiniteRotation(
stage_pole, np.radians(0.01)
)
rotate_slightly_off_mor_opposite_way = (
rotate_slightly_off_mor_one_way.get_inverse()
)
subsegment_index = []
# Iterate over the shared sub-segments.
for (
shared_sub_segment
) in shared_boundary_section.get_shared_sub_segments():
# Tessellate MOR section.
mor_points = pygplates.MultiPointOnSphere(
shared_sub_segment.get_resolved_geometry().to_tessellated(
np.radians(self.ridge_sampling)
)
)
coords = mor_points.to_lat_lon_list()
lats = [i[0] for i in coords]
lons = [i[1] for i in coords]
left_plate = (
shared_boundary_section.get_feature().get_left_plate(None)
)
right_plate = (
shared_boundary_section.get_feature().get_right_plate(None)
)
if left_plate is not None and right_plate is not None:
# Get the spreading rates for all points in this sub segment
(
spreading_rates,
subsegment_index,
) = tools.calculate_spreading_rates(
time=time,
lons=lons,
lats=lats,
left_plates=[left_plate] * len(lons),
right_plates=[right_plate] * len(lons),
rotation_model=self.rotation_model,
delta_time=self.ridge_time_step,
)
else:
spreading_rates = [np.nan] * len(lons)
# Loop through all but the 1st and last points in the current sub segment
for point, rate in zip(
mor_points.get_points()[1:-1],
spreading_rates[1:-1],
):
# Add the point "twice" to the main shifted_mor_points list; once for a L-side
# spread, another for a R-side spread. Then add the same spreading rate twice
# to the list - this therefore assumes spreading rate is symmetric.
shifted_mor_points.append(
rotate_slightly_off_mor_one_way * point
)
shifted_mor_points.append(
rotate_slightly_off_mor_opposite_way * point
)
point_spreading_rates.extend([rate] * 2)
# point_indices.extend(subsegment_index)
# Summarising get_isochrons_for_ridge_snapshot;
# Write out the ridge point born at 'ridge_time' but their position at 'ridge_time - time_step'.
mor_point_features = []
for curr_point in shifted_mor_points:
feature = pygplates.Feature()
feature.set_geometry(curr_point)
feature.set_valid_time(time, -999) # delete - time_step
# feature.set_name(str(spreading_rate))
# feature = set_shapefile_attribute(feature, spreading_rate, "SPREADING_RATE") # make spreading rate a shapefile attribute
mor_point_features.append(feature)
mor_points = pygplates.FeatureCollection(mor_point_features)
# Write MOR points at `time` to gpmlz
basename = "MOR_plus_one_points_{:0.2f}.gpmlz".format(time)
if self.file_collection is not None:
basename = "{}_{}".format(self.file_collection, basename)
mor_points.write(os.path.join(self.save_directory, basename))
# Make sure the max time dataframe is for the initial ocean points only
if time != self._max_time:
self._collect_point_data_in_dataframe(
mor_points, point_spreading_rates, time
)
logger.info(f"Finished building MOR seedpoints at {time} Ma!")
return
def build_all_MOR_seedpoints(self):
"""Resolve mid-ocean ridges for all times between `min_time` and `max_time`, divide them
into points that make up their shared sub-segments. Rotate these points to the left
and right of the ridge using their stage rotation so that they spread from the ridge.
Z-value allocation to each point is done here. In future, a function (like
the spreading rate function) to calculate general z-data will be an input parameter.
Notes
-----
If MOR seed point building is interrupted, progress is safeguarded as long as
`resume_from_checkpoints` is set to `True`.
This assumes that points spread from ridges symmetrically, with the exception of
large ridge jumps at successive timesteps. Therefore, z-values allocated to ridge-emerging
points will appear symmetrical until changes in spreading ridge geometries create
asymmetries.
In future, this will have a checkpoint save feature so that execution
(which occurs during preparation for ReconstructByTopologies and can take several hours)
can be safeguarded against run interruptions.
Parameters
----------
time : float
The time at which to resolve ridge points and stage-rotate them off the ridge.
Returns
-------
mor_point_features : FeatureCollection
All ridge seed points that have emerged from all ridge topologies at `time`.
These points have spread by being slightly rotated away from
ridge locations at `time`.
References
----------
get_mid_ocean_ridge_seedpoints() has been adapted from
https://github.com/siwill22/agegrid-0.1/blob/master/automatic_age_grid_seeding.py#L117.
"""
# If we mustn't overwrite existing files in the `save_directory`, check the status of MOR seeding
# to know where to start/continue seeding
if self.resume_from_checkpoints:
# Check the last MOR seedpoint gpmlz file that was built
checkpointed_MOR_seedpoints = [
s.split("/")[-1]
for s in glob.glob(self.save_directory + "/" + "*MOR_plus_one_points*")
]
try:
# -2 as an index accesses the age (float type), safeguards against identifying numbers in the SeafloorGrid.file_collection string
last_seed_time = np.sort(
[
float(re.findall(r"\d+", s)[-2])
for s in checkpointed_MOR_seedpoints
]
)[0]
# If none were built yet
except:
last_seed_time = "nil"
# If MOR seeding has not started, start it from the top
if last_seed_time == "nil":
time_array = self.time_array
# If the last seed time it could identify is outside the time bounds of the current instance of SeafloorGrid, start
# from the top (this may happen if we use the same save directory for grids for a new set of times)
elif last_seed_time not in self.time_array:
time_array = self.time_array
# If seeding was done to the min_time, we are finished
elif last_seed_time == self.min_time:
return
# If seeding to `min_time` has been interrupted, resume it at last_masked_time.
else:
time_array = np.arange(
last_seed_time, self.min_time - 0.1, -self.ridge_time_step
)
# If we must overwrite all files in `save_directory`, start from `max_time`.
else:
time_array = self.time_array
# Build all continental masks and spreading ridge points (with z values)
self._get_mid_ocean_ridge_seedpoints(time_array)
return
def _create_continental_mask(self, time_array):
"""Create a continental mask for each timestep."""
if time_array[0] != self._max_time:
print(
"Masking interrupted - resuming continental mask building at {} Ma!".format(
time_array[0]
)
)
for time in time_array:
self._PlotTopologies_object.time = time
geoms = self._PlotTopologies_object.continents
final_grid = grids.rasterise(
geoms,
key=1.0,
shape=(self.spacingY, self.spacingX),
extent=self.extent,
origin="lower",
)
final_grid[np.isnan(final_grid)] = 0.0
output_basename = "continent_mask_{}Ma.nc".format(time)
if self.file_collection is not None:
output_basename = "{}_{}".format(
self.file_collection,
output_basename,
)
output_filename = os.path.join(
self.save_directory,
output_basename,
)
grids.write_netcdf_grid(
output_filename, final_grid, extent=[-180, 180, -90, 90]
)
logger.info(f"Finished building a continental mask at {time} Ma!")
return
def build_all_continental_masks(self):
"""Create a continental mask to define the ocean basin for all times between
`min_time` and `max_time`. as well as to use as continental collision
boundaries in `ReconstructByTopologies`.
Notes
-----
Continental masking progress is safeguarded if ever masking is interrupted,
provided that `resume_from_checkpoints` is set to `True`.
If `ReconstructByTopologies` identifies a continental collision
between oceanic points and the boundaries of this continental
mask at `time`, those points are deleted at `time`.
The continental mask is also saved to "/continent_mask_{}Ma.nc" as a
compressed netCDF4 file if a `save_directory` is passed. Otherwise,
the final grid is returned as a NumPy ndarray object.
Returns
-------
all_continental_masks : list of ndarray
A masked grid per timestep in `time_array` with 1=continental point,
and 0=ocean point, for all points on the full global icosahedral mesh.
"""
# If we mustn't overwrite existing files in the `save_directory`, check the status
# of continental masking to know where to start/continue masking
if self.resume_from_checkpoints:
# Check the last continental mask that could be built
checkpointed_continental_masks = [
s.split("/")[-1]
for s in glob.glob(self.save_directory + "/" + "*continent_mask*")
]
try:
# -2 as an index accesses the age (float type), safeguards against identifying numbers in the SeafloorGrid.file_collection string
last_masked_time = np.sort(
[
float(re.findall(r"\d+", s)[-2])
for s in checkpointed_continental_masks
]
)[0]
# If none were built yet
except:
last_masked_time = "nil"
# If masking has not started, start it from the top
if last_masked_time == "nil":
time_array = self.time_array
# If the last seed time it could identify is outside the time bounds of the current instance of SeafloorGrid, start
# from the top (this may happen if we use the same save directory for grids for a new set of times)
elif last_masked_time not in self.time_array:
time_array = self.time_array
# If masking was done to the min_time, we are finished
elif last_masked_time == self.min_time:
return
# If masking to `min_time` has been interrupted, resume it at last_masked_time.
else:
time_array = np.arange(
last_masked_time, self.min_time - 0.1, -self.ridge_time_step
)
# If we must overwrite all files in `save_directory`, start from `max_time`.
else:
time_array = self.time_array
# Build all continental masks and spreading ridge points (with z values)
self._create_continental_mask(time_array)
return
def _extract_zvalues_from_npz_to_ndarray(self, featurecollection, time):
# NPZ file of seedpoint z values that emerged at this time
basename = "point_data_dataframe_{}Ma.npz".format(time)
if self.file_collection is not None:
basename = "{}_{}".format(self.file_collection, basename)
filename = os.path.join(self.save_directory, basename)
loaded_npz = np.load(filename)
curr_zvalues = np.empty([len(featurecollection), len(self.zval_names)])
for i in range(len(self.zval_names)):
# Account for the 0th index being for point feature IDs
curr_zvalues[:, i] = np.array(loaded_npz["arr_{}".format(i)])
return curr_zvalues
def prepare_for_reconstruction_by_topologies(self):
"""Prepare three main auxiliary files for seafloor data gridding:
* Initial ocean seed points (at `max_time`)
* Continental masks (from `max_time` to `min_time`)
* MOR points (from `max_time` to `min_time`)
Returns lists of all attributes for the initial ocean point mesh and
all ridge points for all times in the reconstruction time array.
"""
# INITIAL OCEAN SEED POINT MESH ----------------------------------------------------
(
initial_ocean_seed_points,
initial_ocean_seed_points_mp,
) = self.create_initial_ocean_seed_points()
logger.info("Finished building initial_ocean_seed_points!")
# MOR SEED POINTS AND CONTINENTAL MASKS --------------------------------------------
# The start time for seeding is controlled by whether the overwrite_existing_gridding_inputs
# parameter is set to `True` (in which case the start time is `max_time`). If it is `False`
# and;
# - a run of seeding and continental masking was interrupted, and ridge points were
# checkpointed at n Ma, seeding resumes at n-1 Ma until `min_time` or another interruption
# occurs;
# - seeding was completed but the subsequent gridding input creation was interrupted,
# seeding is assumed completed and skipped. The workflow automatically proceeds to re-gridding.
if self.continent_mask_filename is None:
self.build_all_continental_masks()
else:
logger.info(
"Continent masks passed to SeafloorGrid - skipping continental mask generation!"
)
self.build_all_MOR_seedpoints()
# ALL-TIME POINTS -----------------------------------------------------
# Extract all feature attributes for all reconstruction times into lists
active_points = []
appearance_time = []
birth_lat = [] # latitude_of_crust_formation
prev_lat = []
prev_lon = []
# Extract point feature attributes from MOR seed points
all_mor_features = []
zvalues = np.empty((0, len(self.zval_names)))
for time in self.time_array:
# If we're at the maximum time, start preparing points from the initial ocean mesh
# as well as their z values
if time == self._max_time:
for feature in initial_ocean_seed_points:
active_points.append(feature.get_geometry())
appearance_time.append(feature.get_valid_time()[0])
birth_lat.append(feature.get_geometry().to_lat_lon_list()[0][0])
prev_lat.append(feature.get_geometry().to_lat_lon_list()[0][0])
prev_lon.append(feature.get_geometry().to_lat_lon_list()[0][1])
curr_zvalues = self._extract_zvalues_from_npz_to_ndarray(
initial_ocean_seed_points, time
)
zvalues = np.concatenate((zvalues, curr_zvalues), axis=0)
# Otherwise, we'd be preparing MOR points and their z values
else:
# GPMLZ file of MOR seedpoints
basename = "MOR_plus_one_points_{:0.2f}.gpmlz".format(time)
if self.file_collection is not None:
basename = "{}_{}".format(self.file_collection, basename)
filename = os.path.join(self.save_directory, basename)
features = pygplates.FeatureCollection(filename)
for feature in features:
if feature.get_valid_time()[0] < self.time_array[0]:
active_points.append(feature.get_geometry())
appearance_time.append(feature.get_valid_time()[0])
birth_lat.append(feature.get_geometry().to_lat_lon_list()[0][0])
prev_lat.append(feature.get_geometry().to_lat_lon_list()[0][0])
prev_lon.append(feature.get_geometry().to_lat_lon_list()[0][1])
# COLLECT NDARRAY OF ALL ZVALUES IN THIS TIMESTEP ------------------
curr_zvalues = self._extract_zvalues_from_npz_to_ndarray(features, time)
zvalues = np.concatenate((zvalues, curr_zvalues), axis=0)
return active_points, appearance_time, birth_lat, prev_lat, prev_lon, zvalues
def reconstruct_by_topologies(self):
"""Obtain all active ocean seed points at `time` - these are
points that have not been consumed at subduction zones or have not
collided with continental polygons.
All active points' latitudes, longitues, seafloor ages, spreading rates and all
other general z-values are saved to a gridding input file (.npz).
"""
logger.info("Preparing all initial files...")
# Obtain all info from the ocean seed points and all MOR points through time, store in
# arrays
(
active_points,
appearance_time,
birth_lat,
prev_lat,
prev_lon,
zvalues,
) = self.prepare_for_reconstruction_by_topologies()
#### Begin reconstruction by topology process:
# Indices for all points (`active_points`) that have existed from `max_time` to `min_time`.
point_id = range(len(active_points))
# Specify the default collision detection region as subduction zones
default_collision = reconstruction._DefaultCollision(
feature_specific_collision_parameters=[
(
pygplates.FeatureType.gpml_subduction_zone,
self.subduction_collision_parameters,
)
]
)
# In addition to the default subduction detection, also detect continental collisions
# Use the input continent mask if it is provided.
if self.continent_mask_filename is not None:
collision_spec = reconstruction._ContinentCollision(
# This filename string should not have a time formatted into it - this is
# taken care of later.
self.continent_mask_filename,
default_collision,
verbose=False,
)
else:
# If a continent mask is not provided, use the ones made.
mask_basename = r"continent_mask_{}Ma.nc"
if self.file_collection is not None:
mask_basename = str(self.file_collection) + "_" + mask_basename
mask_template = os.path.join(self.save_directory, mask_basename)
collision_spec = reconstruction._ContinentCollision(
mask_template,
default_collision,
verbose=False,
)
# Call the reconstruct by topologies object
topology_reconstruction = reconstruction._ReconstructByTopologies(
self.rotation_model,
self.topology_features,
self._max_time,
self.min_time,
self.ridge_time_step,
active_points,
point_begin_times=appearance_time,
detect_collisions=collision_spec,
)
# Initialise the reconstruction.
topology_reconstruction.begin_reconstruction()
# Loop over the reconstruction times until the end of the reconstruction time span, or until
# all points have entered their valid time range *and* either exited their time range or
# have been deactivated (subducted forward in time or consumed by MOR backward in time).
reconstruction_data = []
while True:
logger.info(
f"Reconstruct by topologies: working on time {topology_reconstruction.get_current_time():0.2f} Ma"
)
# NOTE:
# topology_reconstruction.get_active_current_points() and topology_reconstruction.get_all_current_points()
# are different. The former is a subset of the latter, and it represents all points at the timestep that
# have not collided with a continental or subduction boundary. The remainders in the latter are inactive
# (NoneType) points, which represent the collided points.
# We need to access active point data from topology_reconstruction.get_all_current_points() because it has
# the same length as the list of all initial ocean points and MOR seed points that have ever emerged from
# spreading ridge topologies through `max_time` to `min_time`. Therefore, it protects the time and space
# order in which all MOR points through time were seeded by pyGPlates. At any given timestep, not all these
# points will be active, but their indices are retained. Thus, z value allocation, point latitudes and
# longitudes of active points will be correctly indexed if taking it from
# topology_reconstruction.get_all_current_points().
curr_points = topology_reconstruction.get_active_current_points()
curr_points_including_inactive = (
topology_reconstruction.get_all_current_points()
)
# Collect latitudes and longitudes of currently ACTIVE points in the ocean basin
curr_lat_lon_points = [point.to_lat_lon() for point in curr_points]
if curr_lat_lon_points:
# Get the number of active points at this timestep.
num_current_points = len(curr_points)
# ndarray to fill with active point lats, lons and zvalues
# FOR NOW, the number of gridding input columns is 6:
# 0 = longitude
# 1 = latitude
# 2 = seafloor age
# 3 = birth latitude snapshot
# 4 = point id
# 5 for the default gridding columns above, plus additional zvalues added next
total_number_of_columns = 5 + len(self.zval_names)
gridding_input_data = np.empty(
[num_current_points, total_number_of_columns]
)
# Lons and lats are first and second columns of the ndarray respectively
gridding_input_data[:, 1], gridding_input_data[:, 0] = zip(
*curr_lat_lon_points
)
# NOTE: We need a single index to access data from curr_points_including_inactive AND allocate
# this data to an ndarray with a number of rows equal to num_current_points. This index will
# append +1 after each loop through curr_points_including_inactive.
i = 0
# Get indices and points of all points at `time`, both active and inactive (which are NoneType points that
# have undergone continental collision or subduction at `time`).
for point_index, current_point in enumerate(
curr_points_including_inactive
):
# Look at all active points (these have not collided with a continent or trench)
if current_point is not None:
# Seafloor age
gridding_input_data[i, 2] = (
appearance_time[point_index]
- topology_reconstruction.get_current_time()
)
# Birth latitude (snapshot)
gridding_input_data[i, 3] = birth_lat[point_index]
# Point ID (snapshot)
gridding_input_data[i, 4] = point_id[
point_index
] # The ID of a corresponding point from the original list of all MOR-resolved points
# GENERAL Z-VALUE ALLOCATION
# Z values are 1st index onwards; 0th belongs to the point feature ID (thus +1)
for j in range(len(self.zval_names)):
# Adjusted index - and we have to add j to 5 to account for lat, lon, age, birth lat and point ID,
adjusted_index = 5 + j
# Spreading rate would be first
# Access current zval from the master list of all zvalues for all points that ever existed in time_array
gridding_input_data[i, adjusted_index] = zvalues[
point_index, j
]
# Go to the next active point
i += 1
gridding_input_dictionary = {}
for i in list(range(total_number_of_columns)):
gridding_input_dictionary[self.total_column_headers[i]] = [
list(j) for j in zip(*gridding_input_data)
][i]
data_to_store = [
gridding_input_dictionary[i] for i in gridding_input_dictionary
]
gridding_input_basename = "gridding_input_{:0.1f}Ma".format(
topology_reconstruction.get_current_time()
)
if self.file_collection is not None:
gridding_input_basename = "{}_{}".format(
self.file_collection,
gridding_input_basename,
)
gridding_input_filename = os.path.join(
self.save_directory, gridding_input_basename
)
# save debug file
if (
"GPLATELY_DEBUG" in os.environ
and os.environ["GPLATELY_DEBUG"].lower() == "true"
):
save_age_grid_sample_points_to_gpml(
gridding_input_dictionary["CURRENT_LONGITUDES"],
gridding_input_dictionary["CURRENT_LATITUDES"],
gridding_input_dictionary["SEAFLOOR_AGE"],
topology_reconstruction.get_current_time(),
self.save_directory,
)
np.savez_compressed(gridding_input_filename, *data_to_store)
if not topology_reconstruction.reconstruct_to_next_time():
break
logger.info(
f"Reconstruction done for {topology_reconstruction.get_current_time()}!"
)
# return reconstruction_data
def lat_lon_z_to_netCDF(
self,
zval_name,
time_arr=None,
unmasked=False,
nprocs=1,
):
"""Produce a netCDF4 grid of a z-value identified by its `zval_name` for a
given time range in `time_arr`.
Seafloor age can be gridded by passing `zval_name` as `SEAFLOOR_AGE`, and spreading
rate can be gridded with `SPREADING_RATE`.
Saves all grids to compressed netCDF format in the attributed directory. Grids
can be read into ndarray format using `gplately.grids.read_netcdf_grid()`.
Parameters
----------
zval_name : str
A string identifiers for a column in the ReconstructByTopologies gridding
input files.
time_arr : list of float, default None
A time range to turn lons, lats and z-values into netCDF4 grids. If not provided,
`time_arr` defaults to the full `time_array` provided to `SeafloorGrids`.
unmasked : bool, default False
Save unmasked grids, in addition to masked versions.
nprocs : int, defaullt 1
Number of processes to use for certain operations (requires joblib).
Passed to `joblib.Parallel`, so -1 means all available processes.
"""
parallel = None
nprocs = int(nprocs)
if nprocs != 1:
try:
from joblib import Parallel
parallel = Parallel(nprocs)
except ImportError:
warnings.warn(
"Could not import joblib; falling back to serial execution"
)
# User can put any time array within SeafloorGrid bounds, but if none
# is provided, it defaults to the attributed time array
if time_arr is None:
time_arr = self.time_array
if parallel is None:
for time in time_arr:
_lat_lon_z_to_netCDF_time(
time=time,
zval_name=zval_name,
file_collection=self.file_collection,
save_directory=self.save_directory,
total_column_headers=self.total_column_headers,
extent=self.extent,
resX=self.spacingX,
resY=self.spacingY,
unmasked=unmasked,
continent_mask_filename=self.continent_mask_filename,
)
else:
from joblib import delayed
parallel(
delayed(_lat_lon_z_to_netCDF_time)(
time=time,
zval_name=zval_name,
file_collection=self.file_collection,
save_directory=self.save_directory,
total_column_headers=self.total_column_headers,
extent=self.extent,
resX=self.spacingX,
resY=self.spacingY,
unmasked=unmasked,
continent_mask_filename=self.continent_mask_filename,
)
for time in time_arr
)
def _lat_lon_z_to_netCDF_time(
time,
zval_name,
file_collection,
save_directory,
total_column_headers,
extent,
resX,
resY,
unmasked=False,
continent_mask_filename=None,
):
# Read the gridding input made by ReconstructByTopologies:
gridding_basename = "gridding_input_{:0.1f}Ma.npz".format(time)
if file_collection is not None:
gridding_basename = "{}_{}".format(file_collection, gridding_basename)
gridding_input = os.path.join(save_directory, gridding_basename)
# Use pandas to load in lons, lats and z values from npz files
npz = np.load(gridding_input)
curr_data = pd.DataFrame.from_dict(
{item: npz[item] for item in npz.files}, orient="columns"
)
curr_data.columns = total_column_headers
# Drop duplicate latitudes and longitudes
unique_data = curr_data.drop_duplicates(
subset=["CURRENT_LONGITUDES", "CURRENT_LATITUDES"]
)
# Acquire lons, lats and zvalues for each time
lons = unique_data["CURRENT_LONGITUDES"].to_list()
lats = unique_data["CURRENT_LATITUDES"].to_list()
zdata = np.array(unique_data[zval_name].to_list())
# zdata = np.where(zdata > 375, float("nan"), zdata), to deal with vmax in the future
zdata = np.nan_to_num(zdata)
# Create a regular grid on which to interpolate lats, lons and zdata
extent_globe = extent
grid_lon = np.linspace(extent_globe[0], extent_globe[1], resX)
grid_lat = np.linspace(extent_globe[2], extent_globe[3], resY)
X, Y = np.meshgrid(grid_lon, grid_lat)
# Interpolate lons, lats and zvals over a regular grid using nearest
# neighbour interpolation
Z = tools.griddata_sphere((lons, lats), zdata, (X, Y), method="nearest")
# Access continental grids from the save directory
if continent_mask_filename is None:
mask_basename = "continent_mask_{}Ma.nc"
if file_collection is not None:
mask_basename = "{}_{}".format(file_collection, mask_basename)
continent_mask_filename = os.path.join(save_directory, mask_basename)
continent_mask_filename = continent_mask_filename.format(time)
unmasked_basename = "{}_grid_unmasked_{}Ma.nc".format(zval_name, time)
grid_basename = "{}_grid_{}Ma.nc".format(zval_name, time)
if file_collection is not None:
unmasked_basename = "{}_{}".format(file_collection, unmasked_basename)
grid_basename = "{}_{}".format(file_collection, grid_basename)
grid_output_unmasked = os.path.join(save_directory, unmasked_basename)
grid_output = os.path.join(save_directory, grid_basename)
if unmasked:
grids.write_netcdf_grid(grid_output_unmasked, Z, extent=extent)
# Identify regions in the grid in the continental mask
cont_mask = grids.Raster(data=continent_mask_filename)
# We need the continental mask to match the number of nodes
# in the uniform grid defined above. This is important if we
# pass our own continental mask to SeafloorGrid
if cont_mask.shape[1] != resX:
cont_mask.resize(resX, resY, inplace=True)
# Use the continental mask
Z = np.ma.array(
grids.Raster(data=Z).data.data, mask=cont_mask.data.data, fill_value=np.nan
)
# grd = cont_mask.interpolate(X, Y) > 0.5
# Z[grd] = np.nan
grids.write_netcdf_grid(
grid_output,
Z,
extent=extent,
)
logger.info(f"netCDF grids for {time} Ma complete!")
def save_age_grid_sample_points_to_gpml(lons, lats, seafloor_ages, paleo_time, outdir):
logger.debug(f"saving age grid sample points to gpml file -- {paleo_time} Ma")
features = []
for lon, lat, age in zip(lons, lats, seafloor_ages):
f = pygplates.Feature()
p = pygplates.PointOnSphere(lat, lon)
f.set_geometry(p)
f.set_valid_time(age + paleo_time, paleo_time)
features.append(f)
pygplates.FeatureCollection(features).write(
os.path.join(outdir, f"age_grid_sample_points_{paleo_time}_Ma.gpmlz")
)Functions
def create_icosahedral_mesh(refinement_levels)-
Define a global point mesh with Stripy's icosahedral triangulated mesh and turn all mesh domains into pyGPlates MultiPointOnSphere types.
This global mesh will be masked with a set of continental or COB terrane polygons to define the ocean basin at a given reconstruction time. The
refinement_levelsinteger is proportional to the resolution of the mesh and the ocean/continent boundary.Parameters
refinement_levels:int- Refine the number of points in the triangulation. The larger the refinement level, the sharper the ocean basin resolution.
Returns
multi_point:instanceof<pygplates.MultiPointOnSphere>- The longitues and latitudes that make up the icosahedral ocean mesh collated into a MultiPointOnSphere object.
icosahedral_global_mesh:instanceof<stripy.spherical_meshes.icosahedral_mesh>- The original global icosahedral triangulated mesh.
Expand source code
def create_icosahedral_mesh(refinement_levels): """Define a global point mesh with Stripy's [icosahedral triangulated mesh](https://github.com/underworldcode/stripy/blob/294354c00dd72e085a018e69c345d9353c6fafef/stripy/spherical_meshes.py#L27) and turn all mesh domains into pyGPlates MultiPointOnSphere types. This global mesh will be masked with a set of continental or COB terrane polygons to define the ocean basin at a given reconstruction time. The `refinement_levels` integer is proportional to the resolution of the mesh and the ocean/continent boundary. Parameters ---------- refinement_levels : int Refine the number of points in the triangulation. The larger the refinement level, the sharper the ocean basin resolution. Returns ------- multi_point : instance of <pygplates.MultiPointOnSphere> The longitues and latitudes that make up the icosahedral ocean mesh collated into a MultiPointOnSphere object. icosahedral_global_mesh : instance of <stripy.spherical_meshes.icosahedral_mesh> The original global icosahedral triangulated mesh. """ import stripy # Create the ocean basin mesh using Stripy's icosahedral spherical mesh icosahedral_global_mesh = stripy.spherical_meshes.icosahedral_mesh( refinement_levels, include_face_points=False, trisection=False, tree=False ) # Get lons and lats of mesh, and turn them into a MultiPointOnSphere lats_arr = np.rad2deg(icosahedral_global_mesh.lats) lons_arr = np.rad2deg(icosahedral_global_mesh.lons) multi_point = pygplates.MultiPointOnSphere(zip(lats_arr, lons_arr)) return multi_point, icosahedral_global_mesh def ensure_polygon_geometry(reconstructed_polygons, rotation_model, time)-
Ensure COB terrane/continental polygon geometries are polygons with reconstruction plate IDs and valid times.
Notes
This step must be done so that the initial set of ocean basin points (the Stripy icosahedral mesh) can be partitioned into plates using each reconstruction plate ID for the given plate
model.This allows for an oceanic point-in-continental polygon query for every identified plate ID. See documentation for
point_in_polygon_routine()for more details.ensure_polygon_geometry()works as follows: COB terrane/continental polygons are assumed to have been reconstructed already inreconstructed_polygons(a list of type). The list contents are turned into a to be ascribed a PolygonOnSpheregeometry, a reconstruction plate ID, and a valid time. Once finished, this feature collection is turned back into a list of instanceand returned. This revert must be completed for compatibility with the subsequent point-in-polygon routine.
Parameters
reconstructed_polygons:listofinstance <pygplates.ReconstructedFeatureGeometry>- If used in
SeafloorGrid, these are automatically obtained from thePlotTopologies.continentsattribute (the reconstructed continental polygons at the current reconstruction time). rotation_model:instanceof<pygplates.RotationModel>- A parameter for turning the
back into a list of instance for compatibility with the point-in-polygon routine.
Expand source code
def ensure_polygon_geometry(reconstructed_polygons, rotation_model, time): """Ensure COB terrane/continental polygon geometries are polygons with reconstruction plate IDs and valid times. Notes ----- This step must be done so that the initial set of ocean basin points (the Stripy icosahedral mesh) can be partitioned into plates using each reconstruction plate ID for the given plate `model`. This allows for an oceanic point-in-continental polygon query for every identified plate ID. See documentation for `point_in_polygon_routine` for more details. `ensure_polygon_geometry` works as follows: COB terrane/continental polygons are assumed to have been reconstructed already in `reconstructed_polygons` (a list of type <pygplates.ReconstructedFeatureGeometry>). The list contents are turned into a <pygplates.FeatureCollection> to be ascribed a `PolygonOnSphere` geometry, a reconstruction plate ID, and a valid time. Once finished, this feature collection is turned back into a list of instance <pygplates.ReconstructedFeatureGeometry> and returned. This revert must be completed for compatibility with the subsequent point-in-polygon routine. Parameters ---------- reconstructed_polygons : list of instance <pygplates.ReconstructedFeatureGeometry> If used in `SeafloorGrid`, these are automatically obtained from the `PlotTopologies.continents` attribute (the reconstructed continental polygons at the current reconstruction time). rotation_model : instance of <pygplates.RotationModel> A parameter for turning the <pygplates.FeatureCollection> back into a list of instance <pygplates.ReconstructedFeatureGeometry> for compatibility with the point-in-polygon routine. """ continent_FeatCol = [] # self._PlotTopologies_object.continents for n in reconstructed_polygons: continent_FeatCol.append(n.get_feature()) polygon_feats = pygplates.FeatureCollection(continent_FeatCol) # From GPRM's force_polygon_geometries(); set feature attributes # like valid times and plate IDs to each masking polygon polygons = [] for feature in polygon_feats: for geom in feature.get_all_geometries(): polygon = pygplates.Feature(feature.get_feature_type()) polygon.set_geometry(pygplates.PolygonOnSphere(geom)) polygon.set_reconstruction_plate_id(feature.get_reconstruction_plate_id()) # Avoid features in COBTerranes with invalid time if feature.get_valid_time()[0] >= feature.get_valid_time()[1]: polygon.set_valid_time( feature.get_valid_time()[0], feature.get_valid_time()[1] ) polygons.append(polygon) cobter_polygon_features = pygplates.FeatureCollection(polygons) # Turn the feature collection back into ReconstructedFeatureGeometry # objects otherwise it will not work with PIP reconstructed_cobter_polygons = [] pygplates.reconstruct( cobter_polygon_features, rotation_model, reconstructed_cobter_polygons, time ) return reconstructed_cobter_polygons def point_in_polygon_routine(multi_point, COB_polygons)-
Perform Plate Tectonic Tools' point in polygon routine to partition points in a
multi_pointMultiPointOnSphere feature based on whether they are inside or outside the polygons inCOB_polygons.Notes
Assuming the
COB_polygonshave passed throughensure_polygon_geometry(), each polygon should have a plate ID assigned to it.This PIP routine serves two purposes for
SeafloorGrid:1) It identifies continental regions in the icosahedral global mesh MultiPointOnSphere feature and 'erases' in-continent oceanic points for the construction of a continental mask at each timestep;
2) It identifies oceanic points in the icosahedral global mesh. These points will be passed to a function that calculates each point's proximity to its nearest MOR segment (if any) within the polygonal domain of its allocated plate ID. Each distance is divided by half the
initial_ocean_mean_spreading_rate(an attribute ofSeafloorGrids) to determine a simplified seafloor age for each point.Number 2) only happens once at the start of the gridding process to momentarily fill the gridding region with initial ocean points that have set ages (albeit not from a plate model file). After multiple time steps of reconstruction, the ocean basin will be filled with new points (with plate-model prescribed ages) that emerge from ridge topologies.
Returns
pygplates.MultiPointOnSphere(points_in_arr) : instance <pygplates.MultiPointOnSphere>- Point features that are within COB terrane polygons.
pygplates.MultiPointOnSphere(points_out_arr) : instance <pygplates.MultiPointOnSphere>- Point features that are outside COB terrane polygons.
zvals:list- A binary list. If an entry is == 0, its corresponing point in the MultiPointOnSphere object is on the ocean. If == 1, the point is in the COB terrane polygon.
Expand source code
def point_in_polygon_routine(multi_point, COB_polygons): """Perform Plate Tectonic Tools' point in polygon routine to partition points in a `multi_point` MultiPointOnSphere feature based on whether they are inside or outside the polygons in `COB_polygons`. Notes ----- Assuming the `COB_polygons` have passed through `ensure_polygon_geometry`, each polygon should have a plate ID assigned to it. This PIP routine serves two purposes for `SeafloorGrid`: 1) It identifies continental regions in the icosahedral global mesh MultiPointOnSphere feature and 'erases' in-continent oceanic points for the construction of a continental mask at each timestep; 2) It identifies oceanic points in the icosahedral global mesh. These points will be passed to a function that calculates each point's proximity to its nearest MOR segment (if any) within the polygonal domain of its allocated plate ID. Each distance is divided by half the `initial_ocean_mean_spreading_rate` (an attribute of `SeafloorGrids`) to determine a simplified seafloor age for each point. Number 2) only happens once at the start of the gridding process to momentarily fill the gridding region with initial ocean points that have set ages (albeit not from a plate model file). After multiple time steps of reconstruction, the ocean basin will be filled with new points (with plate-model prescribed ages) that emerge from ridge topologies. Returns ------- pygplates.MultiPointOnSphere(points_in_arr) : instance <pygplates.MultiPointOnSphere> Point features that are within COB terrane polygons. pygplates.MultiPointOnSphere(points_out_arr) : instance <pygplates.MultiPointOnSphere> Point features that are outside COB terrane polygons. zvals : list A binary list. If an entry is == 0, its corresponing point in the MultiPointOnSphere object is on the ocean. If == 1, the point is in the COB terrane polygon. """ # Convert MultiPointOnSphere to array of PointOnSphere multi_point = np.array(multi_point.get_points(), dtype="object") # Collect reconstructed geometries of continental polygons polygons = np.empty(len(COB_polygons), dtype="object") for ind, i in enumerate(COB_polygons): if isinstance(i, pygplates.ReconstructedFeatureGeometry): geom = i.get_reconstructed_geometry() elif isinstance(i, pygplates.GeometryOnSphere): geom = i else: # e.g. ndarray of coordinates geom = pygplates.PolygonOnSphere(i) polygons[ind] = geom proxies = np.ones(polygons.size) pip_result = ptt.utils.points_in_polygons.find_polygons( multi_point, polygons, proxies, all_polygons=False ) # 1 for points in polygons, None for points outside zvals = np.array( pip_result, dtype="float", ).ravel() zvals[np.isnan(zvals)] = 0.0 zvals = zvals.astype("int") points_in_arr = multi_point[zvals == 1] points_out_arr = multi_point[zvals != 1] return ( pygplates.MultiPointOnSphere(points_in_arr), pygplates.MultiPointOnSphere(points_out_arr), zvals, ) def save_age_grid_sample_points_to_gpml(lons, lats, seafloor_ages, paleo_time, outdir)-
Expand source code
def save_age_grid_sample_points_to_gpml(lons, lats, seafloor_ages, paleo_time, outdir): logger.debug(f"saving age grid sample points to gpml file -- {paleo_time} Ma") features = [] for lon, lat, age in zip(lons, lats, seafloor_ages): f = pygplates.Feature() p = pygplates.PointOnSphere(lat, lon) f.set_geometry(p) f.set_valid_time(age + paleo_time, paleo_time) features.append(f) pygplates.FeatureCollection(features).write( os.path.join(outdir, f"age_grid_sample_points_{paleo_time}_Ma.gpmlz") )
Classes
class SeafloorGrid (PlateReconstruction_object, PlotTopologies_object, max_time, min_time, ridge_time_step, save_directory='agegrids', file_collection=None, refinement_levels=5, ridge_sampling=0.5, extent=(-180, 180, -90, 90), grid_spacing=None, subduction_collision_parameters=(5.0, 10.0), initial_ocean_mean_spreading_rate=75.0, resume_from_checkpoints=False, zval_names=('SPREADING_RATE',), continent_mask_filename=None)-
A class to generate grids that track data atop global ocean basin points (which emerge from mid ocean ridges) through geological time.
Parameters
PlateReconstruction_object:instanceof<gplately.PlateReconstruction>- A GPlately PlateReconstruction object with a
and a containing topology features. PlotTopologies_object:instanceof<gplately.PlotTopologies>- A GPlately PlotTopologies object with a continental polygon or COB terrane polygon file to mask grids with.
max_time:float- The maximum time for age gridding.
min_time:float- The minimum time for age gridding.
ridge_time_step:float- The delta time for resolving ridges (and thus age gridding).
save_directory:str, defaultNone'- The top-level directory to save all outputs to.
file_collection:str, defaultNone- A string to identify the plate model used (will be automated later).
refinement_levels:int, default5- Control the number of points in the icosahedral mesh (higher integer means higher resolution of continent masks).
ridge_sampling:float, default0.5- Spatial resolution (in degrees) at which points that emerge from ridges are tessellated.
extent:listoffloatorint, default[-180.,180.,-90.,90.]- A list containing the mininum longitude, maximum longitude, minimum latitude and maximum latitude extents for all masking and final grids.
grid_spacing:float, defaultNone- The degree spacing/interval with which to space grid points across all masking and
final grids. If
grid_spacingis provided, all grids will use it. If not,grid_spacingdefaults to 0.1. subduction_collision_parameters:len-2 tupleoffloat, default(5.0, 10.0)- A 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my)
- initial_ocean_mean_spreading_rate : float, default 75.
- A spreading rate to uniformly allocate to points that define the initial ocean
- basin. These points will have inaccurate ages, but most of them will be phased
- out after points with plate-model prescribed ages emerge from ridges and spread
- to push them towards collision boundaries (where they are deleted).
resume_from_checkpoints:bool, defaultFalse- If set to
True, and the gridding preparation stage (continental masking and/or ridge seed building) is interrupted, SeafloorGrids will resume gridding preparation from the last successful preparation time. If set toFalse, SeafloorGrids will automatically overwrite all files insave_directoryif re-run after interruption, or normally re-run, thus beginning gridding preparation from scratch.Falsewill be useful if data allocated to the MOR seed points need to be augmented. zval_names:listofstr- A list containing string labels for the z values to attribute to points. Will be used as column headers for z value point dataframes.
continent_mask_filename:str- An optional parameter pointing to the full path to a continental mask for each timestep. Assuming the time is in the filename, i.e. "/path/to/continent_mask_0Ma.nc", it should be passed as "/path/to/continent_mask_{}Ma.nc" with curly brackets. Include decimal formatting if needed.
Expand source code
class SeafloorGrid(object): """A class to generate grids that track data atop global ocean basin points (which emerge from mid ocean ridges) through geological time. Parameters ---------- PlateReconstruction_object : instance of <gplately.PlateReconstruction> A GPlately PlateReconstruction object with a <pygplates.RotationModel> and a <pygplates.FeatureCollection> containing topology features. PlotTopologies_object : instance of <gplately.PlotTopologies> A GPlately PlotTopologies object with a continental polygon or COB terrane polygon file to mask grids with. max_time : float The maximum time for age gridding. min_time : float The minimum time for age gridding. ridge_time_step : float The delta time for resolving ridges (and thus age gridding). save_directory : str, default None' The top-level directory to save all outputs to. file_collection : str, default None A string to identify the plate model used (will be automated later). refinement_levels : int, default 5 Control the number of points in the icosahedral mesh (higher integer means higher resolution of continent masks). ridge_sampling : float, default 0.5 Spatial resolution (in degrees) at which points that emerge from ridges are tessellated. extent : list of float or int, default [-180.,180.,-90.,90.] A list containing the mininum longitude, maximum longitude, minimum latitude and maximum latitude extents for all masking and final grids. grid_spacing : float, default None The degree spacing/interval with which to space grid points across all masking and final grids. If `grid_spacing` is provided, all grids will use it. If not, `grid_spacing` defaults to 0.1. subduction_collision_parameters : len-2 tuple of float, default (5.0, 10.0) A 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my) initial_ocean_mean_spreading_rate : float, default 75. A spreading rate to uniformly allocate to points that define the initial ocean basin. These points will have inaccurate ages, but most of them will be phased out after points with plate-model prescribed ages emerge from ridges and spread to push them towards collision boundaries (where they are deleted). resume_from_checkpoints : bool, default False If set to `True`, and the gridding preparation stage (continental masking and/or ridge seed building) is interrupted, SeafloorGrids will resume gridding preparation from the last successful preparation time. If set to `False`, SeafloorGrids will automatically overwrite all files in `save_directory` if re-run after interruption, or normally re-run, thus beginning gridding preparation from scratch. `False` will be useful if data allocated to the MOR seed points need to be augmented. zval_names : list of str A list containing string labels for the z values to attribute to points. Will be used as column headers for z value point dataframes. continent_mask_filename : str An optional parameter pointing to the full path to a continental mask for each timestep. Assuming the time is in the filename, i.e. "/path/to/continent_mask_0Ma.nc", it should be passed as "/path/to/continent_mask_{}Ma.nc" with curly brackets. Include decimal formatting if needed. """ def __init__( self, PlateReconstruction_object, PlotTopologies_object, max_time, min_time, ridge_time_step, save_directory="agegrids", file_collection=None, refinement_levels=5, ridge_sampling=0.5, extent=(-180, 180, -90, 90), grid_spacing=None, subduction_collision_parameters=(5.0, 10.0), initial_ocean_mean_spreading_rate=75.0, resume_from_checkpoints=False, zval_names=("SPREADING_RATE",), continent_mask_filename=None, ): # Provides a rotation model, topology features and reconstruction time for # the SeafloorGrid self.PlateReconstruction_object = PlateReconstruction_object self.rotation_model = self.PlateReconstruction_object.rotation_model self.topology_features = self.PlateReconstruction_object.topology_features self._PlotTopologies_object = PlotTopologies_object save_directory = str(save_directory) if not os.path.isdir(save_directory): logger.info( "Output directory does not exist; creating now: " + save_directory ) os.makedirs(save_directory, exist_ok=True) self.save_directory = save_directory if file_collection is not None: file_collection = str(file_collection) self.file_collection = file_collection # Topological parameters self.refinement_levels = refinement_levels self.ridge_sampling = ridge_sampling self.subduction_collision_parameters = subduction_collision_parameters self.initial_ocean_mean_spreading_rate = initial_ocean_mean_spreading_rate # Gridding parameters self.extent = extent # A list of degree spacings that allow an even division of the global lat-lon extent. divisible_degree_spacings = [0.1, 0.25, 0.5, 0.75, 1.0] if grid_spacing: # If the provided degree spacing is in the list of permissible spacings, use it # and prepare the number of pixels in x and y (spacingX and spacingY) if grid_spacing in divisible_degree_spacings: self.grid_spacing = grid_spacing self.spacingX = _deg2pixels( grid_spacing, self.extent[0], self.extent[1] ) self.spacingY = _deg2pixels( grid_spacing, self.extent[2], self.extent[3] ) # If the provided spacing is >>1 degree, use 1 degree elif grid_spacing >= divisible_degree_spacings[-1]: self.grid_spacing = divisible_degree_spacings[-1] self.spacingX = _deg2pixels( divisible_degree_spacings[-1], self.extent[0], self.extent[1] ) self.spacingY = _deg2pixels( divisible_degree_spacings[-1], self.extent[2], self.extent[3] ) with warnings.catch_warnings(): warnings.simplefilter("always") warnings.warn( f"The provided grid_spacing of {grid_spacing} is quite large. To preserve the grid resolution, a {self.grid_spacing} degree spacing has been employed instead." ) # If the provided degree spacing is not in the list of permissible spacings, but below # a degree, find the closest permissible degree spacing. Use this and find # spacingX and spacingY. else: for divisible_degree_spacing in divisible_degree_spacings: # The tolerance is half the difference between consecutive divisible spacings. # Max is 1 degree for now - other integers work but may provide too little of a # grid resolution. if abs(grid_spacing - divisible_degree_spacing) <= 0.125: new_deg_res = divisible_degree_spacing self.grid_spacing = new_deg_res self.spacingX = _deg2pixels( new_deg_res, self.extent[0], self.extent[1] ) self.spacingY = _deg2pixels( new_deg_res, self.extent[2], self.extent[3] ) with warnings.catch_warnings(): warnings.simplefilter("always") warnings.warn( f"The provided grid_spacing of {grid_spacing} does not cleanly divide into the global extent. A degree spacing of {self.grid_spacing} has been employed instead." ) else: # If a spacing degree is not provided, use default # resolution and get default spacingX and spacingY self.grid_spacing = 0.1 self.spacingX = 3601 self.spacingY = 1801 self.resume_from_checkpoints = resume_from_checkpoints # Temporal parameters self._max_time = float(max_time) self.min_time = float(min_time) self.ridge_time_step = float(ridge_time_step) self.time_array = np.arange( self._max_time, self.min_time - 0.1, -self.ridge_time_step ) # If PlotTopologies' time attribute is not equal to the maximum time in the # seafloor grid reconstruction tree, make it equal. This will ensure the time # for continental masking is consistent. if self._PlotTopologies_object.time != self._max_time: self._PlotTopologies_object.time = self._max_time # Essential features and meshes for the SeafloorGrid self.continental_polygons = ensure_polygon_geometry( self._PlotTopologies_object.continents, self.rotation_model, self._max_time ) self._PlotTopologies_object.continents = PlotTopologies_object.continents ( self.icosahedral_multi_point, self.icosahedral_global_mesh, ) = create_icosahedral_mesh(self.refinement_levels) # Z value parameters self.zval_names = zval_names self.default_column_headers = [ "CURRENT_LONGITUDES", "CURRENT_LATITUDES", "SEAFLOOR_AGE", "BIRTH_LAT_SNAPSHOT", "POINT_ID_SNAPSHOT", ] self.total_column_headers = np.concatenate( [self.default_column_headers, self.zval_names] ) # Filename for continental masks that the user can provide instead of building it here self.continent_mask_filename = continent_mask_filename # If the user provides a continental mask filename, we need to downsize the mask # resolution for when we create the initial ocean mesh. The mesh does not need to be high-res. if self.continent_mask_filename is not None: # Determine which percentage to use to scale the continent mask resolution at max time def _map_res_to_node_percentage(self, continent_mask_filename): maskY, maskX = grids.read_netcdf_grid( continent_mask_filename.format(self._max_time) ).shape mask_deg = _pixels2deg(maskX, self.extent[0], self.extent[1]) if mask_deg <= 0.1: percentage = 0.1 elif mask_deg <= 0.25: percentage = 0.3 elif mask_deg <= 0.5: percentage = 0.5 elif mask_deg < 0.75: percentage = 0.6 elif mask_deg >= 1: percentage = 0.75 return mask_deg, percentage _, self.percentage = _map_res_to_node_percentage( self, self.continent_mask_filename ) # Allow SeafloorGrid time to be updated, and to update the internally-used # PlotTopologies' time attribute too. If PlotTopologies is used outside the # object, its `time` attribute is not updated. @property def max_time(self): """The reconstruction time.""" return self._max_time @property def PlotTopologiesTime(self): return self._PlotTopologies_object.time @max_time.setter def max_time(self, var): if var >= 0: self.update_time(var) else: raise ValueError("Enter a valid time >= 0") def update_time(self, max_time): self._max_time = float(max_time) self._PlotTopologies_object.time = float(max_time) def _collect_point_data_in_dataframe( self, pygplates_featurecollection, zval_ndarray, time ): """At a given timestep, create a pandas dataframe holding all attributes of point features. Rather than store z values as shapefile attributes, store them in a dataframe indexed by feature ID. """ # Turn the zval_ndarray into a numPy array zval_ndarray = np.array(zval_ndarray) feature_id = [] for feature in pygplates_featurecollection: feature_id.append(str(feature.get_feature_id())) # Prepare the zval ndarray (can be of any shape) to be saved with default point data zvals_to_store = {} # If only one zvalue (fow now, spreading rate) if zval_ndarray.ndim == 1: zvals_to_store[self.zval_names[0]] = zval_ndarray data_to_store = [zvals_to_store[i] for i in zvals_to_store] else: for i in zval_ndarray.shape[1]: zvals_to_store[self.zval_names[i]] = [ list(j) for j in zip(*zval_ndarray) ][i] data_to_store = [zvals_to_store[i] for i in zvals_to_store] basename = "point_data_dataframe_{}Ma".format(time) if self.file_collection is not None: basename = "{}_{}".format(self.file_collection, basename) filename = os.path.join(self.save_directory, basename) np.savez_compressed(filename, FEATURE_ID=feature_id, *data_to_store) return def create_initial_ocean_seed_points(self): """Create the initial ocean basin seed point domain (at `max_time` only) using Stripy's icosahedral triangulation with the specified `self.refinement_levels`. The ocean mesh starts off as a global-spanning Stripy icosahedral mesh. `create_initial_ocean_seed_points` passes the automatically-resolved-to-current-time continental polygons from the `PlotTopologies_object`'s `continents` attribute (which can be from a COB terrane file or a continental polygon file) into Plate Tectonic Tools' point-in-polygon routine. It identifies ocean basin points that lie: * outside the polygons (for the ocean basin point domain) * inside the polygons (for the continental mask) Points from the mesh outside the continental polygons make up the ocean basin seed point mesh. The masked mesh is outputted as a compressed GPML (GPMLZ) file with the filename: "ocean_basin_seed_points_{}Ma.gpmlz" if a `save_directory` is passed. Otherwise, the mesh is returned as a pyGPlates FeatureCollection object. Notes ----- This point mesh represents ocean basin seafloor that was produced before `SeafloorGrid.max_time`, and thus has unknown properties like valid time and spreading rate. As time passes, the plate reconstruction model sees points emerging from MORs. These new points spread to occupy the ocean basins, moving the initial filler points closer to subduction zones and continental polygons with which they can collide. If a collision is detected by `PlateReconstruction`s `ReconstructByTopologies` object, these points are deleted. Ideally, if a reconstruction tree spans a large time range, **all** initial mesh points would collide with a continent or be subducted, leaving behind a mesh of well-defined MOR-emerged ocean basin points that data can be attributed to. However, some of these initial points situated close to contiental boundaries are retained through time - these form point artefacts with anomalously high ages. Even deep-time plate models (e.g. 1 Ga) will have these artefacts - removing them would require more detail to be added to the reconstruction model. Returns ------- ocean_basin_point_mesh : pygplates.FeatureCollection of pygplates.MultiPointOnSphere A feature collection of point objects on the ocean basin. """ if self.continent_mask_filename is None: # Ensure COB terranes at max time have reconstruction IDs and valid times COB_polygons = ensure_polygon_geometry( self._PlotTopologies_object.continents, self.rotation_model, self._max_time, ) # zval is a binary array encoding whether a point # coordinate is within a COB terrane polygon or not. # Use the icosahedral mesh MultiPointOnSphere attribute _, ocean_basin_point_mesh, zvals = point_in_polygon_routine( self.icosahedral_multi_point, COB_polygons ) # Plates to partition with plate_partitioner = pygplates.PlatePartitioner( COB_polygons, self.rotation_model, ) # Plate partition the ocean basin points meshnode_feature = pygplates.Feature( pygplates.FeatureType.create_from_qualified_string("gpml:MeshNode") ) meshnode_feature.set_geometry( ocean_basin_point_mesh # multi_point ) ocean_basin_meshnode = pygplates.FeatureCollection(meshnode_feature) paleogeography = plate_partitioner.partition_features( ocean_basin_meshnode, partition_return=pygplates.PartitionReturn.separate_partitioned_and_unpartitioned, properties_to_copy=[pygplates.PropertyName.gpml_shapefile_attributes], ) ocean_points = paleogeography[1] # Separate those inside polygons continent_points = paleogeography[0] # Separate those outside polygons # If a set of continent masks was passed, we can use max_time's continental # mask to build the initial profile of seafloor age. else: max_time_cont_mask = grids.Raster( self.continent_mask_filename.format(self._max_time) ) # If the input grid is at 0.5 degree uniform spacing, then the input # grid is 7x more populated than a 6-level stripy icosahedral mesh and # using this resolution for the initial ocean mesh will dramatically slow down # reconstruction by topologies. # Scale down the resolution based on the input mask resolution # (percentage was found in __init__.) max_time_cont_mask.resize( int(max_time_cont_mask.shape[0] * self.percentage), int(max_time_cont_mask.shape[1] * self.percentage), inplace=True, ) lat = np.linspace(-90, 90, max_time_cont_mask.shape[0]) lon = np.linspace(-180, 180, max_time_cont_mask.shape[1]) llon, llat = np.meshgrid(lon, lat) mask_inds = np.where(max_time_cont_mask.data.flatten() == 0) mask_vals = max_time_cont_mask.data.flatten() mask_lon = llon.flatten()[mask_inds] mask_lat = llat.flatten()[mask_inds] ocean_pt_feature = pygplates.Feature() ocean_pt_feature.set_geometry( pygplates.MultiPointOnSphere(zip(mask_lat, mask_lon)) ) ocean_points = [ocean_pt_feature] # Now that we have ocean points... # Determine age of ocean basin points using their proximity to MOR features # and an assumed globally-uniform ocean basin mean spreading rate. # We need resolved topologies at the `max_time` to pass into the proximity # function resolved_topologies = [] shared_boundary_sections = [] pygplates.resolve_topologies( self.topology_features, self.rotation_model, resolved_topologies, self._max_time, shared_boundary_sections, ) pX, pY, pZ = tools.find_distance_to_nearest_ridge( resolved_topologies, shared_boundary_sections, ocean_points, ) # Divide spreading rate by 2 to use half the mean spreading rate pAge = np.array(pZ) / (self.initial_ocean_mean_spreading_rate / 2.0) initial_ocean_point_features = [] initial_ocean_multipoints = [] for point in zip(pX, pY, pAge): point_feature = pygplates.Feature() point_feature.set_geometry(pygplates.PointOnSphere(point[1], point[0])) # Add 'time' to the age at the time of computation, to get the valid time in Ma point_feature.set_valid_time(point[2] + self._max_time, -1) # For now: custom zvals are added as shapefile attributes - will attempt pandas data frames # point_feature = set_shapefile_attribute(point_feature, self.initial_ocean_mean_spreading_rate, "SPREADING_RATE") # Seems like static data initial_ocean_point_features.append(point_feature) initial_ocean_multipoints.append(point_feature.get_geometry()) # print(initial_ocean_point_features) multi_point_feature = pygplates.MultiPointOnSphere(initial_ocean_multipoints) basename = "ocean_basin_seed_points_{}_RLs_{}Ma.gpmlz".format( self.refinement_levels, self._max_time, ) if self.file_collection is not None: basename = "{}_{}".format(self.file_collection, basename) output_filename = os.path.join(self.save_directory, basename) initial_ocean_feature_collection = pygplates.FeatureCollection( initial_ocean_point_features ) initial_ocean_feature_collection.write(output_filename) # Collect all point feature data into a pandas dataframe self._collect_point_data_in_dataframe( initial_ocean_feature_collection, np.array( [self.initial_ocean_mean_spreading_rate] * len(pX) ), # for now, spreading rate is one zvalue for initial ocean points. will other zvalues need to have a generalised workflow? self._max_time, ) return ( pygplates.FeatureCollection(initial_ocean_point_features), multi_point_feature, ) def _get_mid_ocean_ridge_seedpoints(self, time_array): # Topology features from `PlotTopologies`. topology_features_extracted = pygplates.FeaturesFunctionArgument( self.topology_features ) # Create a mask for each timestep if time_array[0] != self._max_time: print( "MOR seed point building interrupted - resuming at {} Ma!".format( time_array[0] ) ) for time in time_array: # Points and their z values that emerge from MORs at this time. shifted_mor_points = [] point_spreading_rates = [] # Resolve topologies to the current time. resolved_topologies = [] shared_boundary_sections = [] pygplates.resolve_topologies( topology_features_extracted.get_features(), self.rotation_model, resolved_topologies, time, shared_boundary_sections, ) # pygplates.ResolvedTopologicalSection objects. for shared_boundary_section in shared_boundary_sections: if ( shared_boundary_section.get_feature().get_feature_type() == pygplates.FeatureType.create_gpml("MidOceanRidge") ): spreading_feature = shared_boundary_section.get_feature() # Find the stage rotation of the spreading feature in the # frame of reference of its geometry at the current # reconstruction time (the MOR is currently actively spreading). # The stage pole can then be directly geometrically compared # to the *reconstructed* spreading geometry. stage_rotation = separate_ridge_transform_segments.get_stage_rotation_for_reconstructed_geometry( spreading_feature, self.rotation_model, time ) if not stage_rotation: # Skip current feature - it's not a spreading feature. continue # Get the stage pole of the stage rotation. # Note that the stage rotation is already in frame of # reference of the *reconstructed* geometry at the spreading time. stage_pole, _ = stage_rotation.get_euler_pole_and_angle() # One way rotates left and the other right, but don't know # which - doesn't matter in our example though. rotate_slightly_off_mor_one_way = pygplates.FiniteRotation( stage_pole, np.radians(0.01) ) rotate_slightly_off_mor_opposite_way = ( rotate_slightly_off_mor_one_way.get_inverse() ) subsegment_index = [] # Iterate over the shared sub-segments. for ( shared_sub_segment ) in shared_boundary_section.get_shared_sub_segments(): # Tessellate MOR section. mor_points = pygplates.MultiPointOnSphere( shared_sub_segment.get_resolved_geometry().to_tessellated( np.radians(self.ridge_sampling) ) ) coords = mor_points.to_lat_lon_list() lats = [i[0] for i in coords] lons = [i[1] for i in coords] left_plate = ( shared_boundary_section.get_feature().get_left_plate(None) ) right_plate = ( shared_boundary_section.get_feature().get_right_plate(None) ) if left_plate is not None and right_plate is not None: # Get the spreading rates for all points in this sub segment ( spreading_rates, subsegment_index, ) = tools.calculate_spreading_rates( time=time, lons=lons, lats=lats, left_plates=[left_plate] * len(lons), right_plates=[right_plate] * len(lons), rotation_model=self.rotation_model, delta_time=self.ridge_time_step, ) else: spreading_rates = [np.nan] * len(lons) # Loop through all but the 1st and last points in the current sub segment for point, rate in zip( mor_points.get_points()[1:-1], spreading_rates[1:-1], ): # Add the point "twice" to the main shifted_mor_points list; once for a L-side # spread, another for a R-side spread. Then add the same spreading rate twice # to the list - this therefore assumes spreading rate is symmetric. shifted_mor_points.append( rotate_slightly_off_mor_one_way * point ) shifted_mor_points.append( rotate_slightly_off_mor_opposite_way * point ) point_spreading_rates.extend([rate] * 2) # point_indices.extend(subsegment_index) # Summarising get_isochrons_for_ridge_snapshot; # Write out the ridge point born at 'ridge_time' but their position at 'ridge_time - time_step'. mor_point_features = [] for curr_point in shifted_mor_points: feature = pygplates.Feature() feature.set_geometry(curr_point) feature.set_valid_time(time, -999) # delete - time_step # feature.set_name(str(spreading_rate)) # feature = set_shapefile_attribute(feature, spreading_rate, "SPREADING_RATE") # make spreading rate a shapefile attribute mor_point_features.append(feature) mor_points = pygplates.FeatureCollection(mor_point_features) # Write MOR points at `time` to gpmlz basename = "MOR_plus_one_points_{:0.2f}.gpmlz".format(time) if self.file_collection is not None: basename = "{}_{}".format(self.file_collection, basename) mor_points.write(os.path.join(self.save_directory, basename)) # Make sure the max time dataframe is for the initial ocean points only if time != self._max_time: self._collect_point_data_in_dataframe( mor_points, point_spreading_rates, time ) logger.info(f"Finished building MOR seedpoints at {time} Ma!") return def build_all_MOR_seedpoints(self): """Resolve mid-ocean ridges for all times between `min_time` and `max_time`, divide them into points that make up their shared sub-segments. Rotate these points to the left and right of the ridge using their stage rotation so that they spread from the ridge. Z-value allocation to each point is done here. In future, a function (like the spreading rate function) to calculate general z-data will be an input parameter. Notes ----- If MOR seed point building is interrupted, progress is safeguarded as long as `resume_from_checkpoints` is set to `True`. This assumes that points spread from ridges symmetrically, with the exception of large ridge jumps at successive timesteps. Therefore, z-values allocated to ridge-emerging points will appear symmetrical until changes in spreading ridge geometries create asymmetries. In future, this will have a checkpoint save feature so that execution (which occurs during preparation for ReconstructByTopologies and can take several hours) can be safeguarded against run interruptions. Parameters ---------- time : float The time at which to resolve ridge points and stage-rotate them off the ridge. Returns ------- mor_point_features : FeatureCollection All ridge seed points that have emerged from all ridge topologies at `time`. These points have spread by being slightly rotated away from ridge locations at `time`. References ---------- get_mid_ocean_ridge_seedpoints() has been adapted from https://github.com/siwill22/agegrid-0.1/blob/master/automatic_age_grid_seeding.py#L117. """ # If we mustn't overwrite existing files in the `save_directory`, check the status of MOR seeding # to know where to start/continue seeding if self.resume_from_checkpoints: # Check the last MOR seedpoint gpmlz file that was built checkpointed_MOR_seedpoints = [ s.split("/")[-1] for s in glob.glob(self.save_directory + "/" + "*MOR_plus_one_points*") ] try: # -2 as an index accesses the age (float type), safeguards against identifying numbers in the SeafloorGrid.file_collection string last_seed_time = np.sort( [ float(re.findall(r"\d+", s)[-2]) for s in checkpointed_MOR_seedpoints ] )[0] # If none were built yet except: last_seed_time = "nil" # If MOR seeding has not started, start it from the top if last_seed_time == "nil": time_array = self.time_array # If the last seed time it could identify is outside the time bounds of the current instance of SeafloorGrid, start # from the top (this may happen if we use the same save directory for grids for a new set of times) elif last_seed_time not in self.time_array: time_array = self.time_array # If seeding was done to the min_time, we are finished elif last_seed_time == self.min_time: return # If seeding to `min_time` has been interrupted, resume it at last_masked_time. else: time_array = np.arange( last_seed_time, self.min_time - 0.1, -self.ridge_time_step ) # If we must overwrite all files in `save_directory`, start from `max_time`. else: time_array = self.time_array # Build all continental masks and spreading ridge points (with z values) self._get_mid_ocean_ridge_seedpoints(time_array) return def _create_continental_mask(self, time_array): """Create a continental mask for each timestep.""" if time_array[0] != self._max_time: print( "Masking interrupted - resuming continental mask building at {} Ma!".format( time_array[0] ) ) for time in time_array: self._PlotTopologies_object.time = time geoms = self._PlotTopologies_object.continents final_grid = grids.rasterise( geoms, key=1.0, shape=(self.spacingY, self.spacingX), extent=self.extent, origin="lower", ) final_grid[np.isnan(final_grid)] = 0.0 output_basename = "continent_mask_{}Ma.nc".format(time) if self.file_collection is not None: output_basename = "{}_{}".format( self.file_collection, output_basename, ) output_filename = os.path.join( self.save_directory, output_basename, ) grids.write_netcdf_grid( output_filename, final_grid, extent=[-180, 180, -90, 90] ) logger.info(f"Finished building a continental mask at {time} Ma!") return def build_all_continental_masks(self): """Create a continental mask to define the ocean basin for all times between `min_time` and `max_time`. as well as to use as continental collision boundaries in `ReconstructByTopologies`. Notes ----- Continental masking progress is safeguarded if ever masking is interrupted, provided that `resume_from_checkpoints` is set to `True`. If `ReconstructByTopologies` identifies a continental collision between oceanic points and the boundaries of this continental mask at `time`, those points are deleted at `time`. The continental mask is also saved to "/continent_mask_{}Ma.nc" as a compressed netCDF4 file if a `save_directory` is passed. Otherwise, the final grid is returned as a NumPy ndarray object. Returns ------- all_continental_masks : list of ndarray A masked grid per timestep in `time_array` with 1=continental point, and 0=ocean point, for all points on the full global icosahedral mesh. """ # If we mustn't overwrite existing files in the `save_directory`, check the status # of continental masking to know where to start/continue masking if self.resume_from_checkpoints: # Check the last continental mask that could be built checkpointed_continental_masks = [ s.split("/")[-1] for s in glob.glob(self.save_directory + "/" + "*continent_mask*") ] try: # -2 as an index accesses the age (float type), safeguards against identifying numbers in the SeafloorGrid.file_collection string last_masked_time = np.sort( [ float(re.findall(r"\d+", s)[-2]) for s in checkpointed_continental_masks ] )[0] # If none were built yet except: last_masked_time = "nil" # If masking has not started, start it from the top if last_masked_time == "nil": time_array = self.time_array # If the last seed time it could identify is outside the time bounds of the current instance of SeafloorGrid, start # from the top (this may happen if we use the same save directory for grids for a new set of times) elif last_masked_time not in self.time_array: time_array = self.time_array # If masking was done to the min_time, we are finished elif last_masked_time == self.min_time: return # If masking to `min_time` has been interrupted, resume it at last_masked_time. else: time_array = np.arange( last_masked_time, self.min_time - 0.1, -self.ridge_time_step ) # If we must overwrite all files in `save_directory`, start from `max_time`. else: time_array = self.time_array # Build all continental masks and spreading ridge points (with z values) self._create_continental_mask(time_array) return def _extract_zvalues_from_npz_to_ndarray(self, featurecollection, time): # NPZ file of seedpoint z values that emerged at this time basename = "point_data_dataframe_{}Ma.npz".format(time) if self.file_collection is not None: basename = "{}_{}".format(self.file_collection, basename) filename = os.path.join(self.save_directory, basename) loaded_npz = np.load(filename) curr_zvalues = np.empty([len(featurecollection), len(self.zval_names)]) for i in range(len(self.zval_names)): # Account for the 0th index being for point feature IDs curr_zvalues[:, i] = np.array(loaded_npz["arr_{}".format(i)]) return curr_zvalues def prepare_for_reconstruction_by_topologies(self): """Prepare three main auxiliary files for seafloor data gridding: * Initial ocean seed points (at `max_time`) * Continental masks (from `max_time` to `min_time`) * MOR points (from `max_time` to `min_time`) Returns lists of all attributes for the initial ocean point mesh and all ridge points for all times in the reconstruction time array. """ # INITIAL OCEAN SEED POINT MESH ---------------------------------------------------- ( initial_ocean_seed_points, initial_ocean_seed_points_mp, ) = self.create_initial_ocean_seed_points() logger.info("Finished building initial_ocean_seed_points!") # MOR SEED POINTS AND CONTINENTAL MASKS -------------------------------------------- # The start time for seeding is controlled by whether the overwrite_existing_gridding_inputs # parameter is set to `True` (in which case the start time is `max_time`). If it is `False` # and; # - a run of seeding and continental masking was interrupted, and ridge points were # checkpointed at n Ma, seeding resumes at n-1 Ma until `min_time` or another interruption # occurs; # - seeding was completed but the subsequent gridding input creation was interrupted, # seeding is assumed completed and skipped. The workflow automatically proceeds to re-gridding. if self.continent_mask_filename is None: self.build_all_continental_masks() else: logger.info( "Continent masks passed to SeafloorGrid - skipping continental mask generation!" ) self.build_all_MOR_seedpoints() # ALL-TIME POINTS ----------------------------------------------------- # Extract all feature attributes for all reconstruction times into lists active_points = [] appearance_time = [] birth_lat = [] # latitude_of_crust_formation prev_lat = [] prev_lon = [] # Extract point feature attributes from MOR seed points all_mor_features = [] zvalues = np.empty((0, len(self.zval_names))) for time in self.time_array: # If we're at the maximum time, start preparing points from the initial ocean mesh # as well as their z values if time == self._max_time: for feature in initial_ocean_seed_points: active_points.append(feature.get_geometry()) appearance_time.append(feature.get_valid_time()[0]) birth_lat.append(feature.get_geometry().to_lat_lon_list()[0][0]) prev_lat.append(feature.get_geometry().to_lat_lon_list()[0][0]) prev_lon.append(feature.get_geometry().to_lat_lon_list()[0][1]) curr_zvalues = self._extract_zvalues_from_npz_to_ndarray( initial_ocean_seed_points, time ) zvalues = np.concatenate((zvalues, curr_zvalues), axis=0) # Otherwise, we'd be preparing MOR points and their z values else: # GPMLZ file of MOR seedpoints basename = "MOR_plus_one_points_{:0.2f}.gpmlz".format(time) if self.file_collection is not None: basename = "{}_{}".format(self.file_collection, basename) filename = os.path.join(self.save_directory, basename) features = pygplates.FeatureCollection(filename) for feature in features: if feature.get_valid_time()[0] < self.time_array[0]: active_points.append(feature.get_geometry()) appearance_time.append(feature.get_valid_time()[0]) birth_lat.append(feature.get_geometry().to_lat_lon_list()[0][0]) prev_lat.append(feature.get_geometry().to_lat_lon_list()[0][0]) prev_lon.append(feature.get_geometry().to_lat_lon_list()[0][1]) # COLLECT NDARRAY OF ALL ZVALUES IN THIS TIMESTEP ------------------ curr_zvalues = self._extract_zvalues_from_npz_to_ndarray(features, time) zvalues = np.concatenate((zvalues, curr_zvalues), axis=0) return active_points, appearance_time, birth_lat, prev_lat, prev_lon, zvalues def reconstruct_by_topologies(self): """Obtain all active ocean seed points at `time` - these are points that have not been consumed at subduction zones or have not collided with continental polygons. All active points' latitudes, longitues, seafloor ages, spreading rates and all other general z-values are saved to a gridding input file (.npz). """ logger.info("Preparing all initial files...") # Obtain all info from the ocean seed points and all MOR points through time, store in # arrays ( active_points, appearance_time, birth_lat, prev_lat, prev_lon, zvalues, ) = self.prepare_for_reconstruction_by_topologies() #### Begin reconstruction by topology process: # Indices for all points (`active_points`) that have existed from `max_time` to `min_time`. point_id = range(len(active_points)) # Specify the default collision detection region as subduction zones default_collision = reconstruction._DefaultCollision( feature_specific_collision_parameters=[ ( pygplates.FeatureType.gpml_subduction_zone, self.subduction_collision_parameters, ) ] ) # In addition to the default subduction detection, also detect continental collisions # Use the input continent mask if it is provided. if self.continent_mask_filename is not None: collision_spec = reconstruction._ContinentCollision( # This filename string should not have a time formatted into it - this is # taken care of later. self.continent_mask_filename, default_collision, verbose=False, ) else: # If a continent mask is not provided, use the ones made. mask_basename = r"continent_mask_{}Ma.nc" if self.file_collection is not None: mask_basename = str(self.file_collection) + "_" + mask_basename mask_template = os.path.join(self.save_directory, mask_basename) collision_spec = reconstruction._ContinentCollision( mask_template, default_collision, verbose=False, ) # Call the reconstruct by topologies object topology_reconstruction = reconstruction._ReconstructByTopologies( self.rotation_model, self.topology_features, self._max_time, self.min_time, self.ridge_time_step, active_points, point_begin_times=appearance_time, detect_collisions=collision_spec, ) # Initialise the reconstruction. topology_reconstruction.begin_reconstruction() # Loop over the reconstruction times until the end of the reconstruction time span, or until # all points have entered their valid time range *and* either exited their time range or # have been deactivated (subducted forward in time or consumed by MOR backward in time). reconstruction_data = [] while True: logger.info( f"Reconstruct by topologies: working on time {topology_reconstruction.get_current_time():0.2f} Ma" ) # NOTE: # topology_reconstruction.get_active_current_points() and topology_reconstruction.get_all_current_points() # are different. The former is a subset of the latter, and it represents all points at the timestep that # have not collided with a continental or subduction boundary. The remainders in the latter are inactive # (NoneType) points, which represent the collided points. # We need to access active point data from topology_reconstruction.get_all_current_points() because it has # the same length as the list of all initial ocean points and MOR seed points that have ever emerged from # spreading ridge topologies through `max_time` to `min_time`. Therefore, it protects the time and space # order in which all MOR points through time were seeded by pyGPlates. At any given timestep, not all these # points will be active, but their indices are retained. Thus, z value allocation, point latitudes and # longitudes of active points will be correctly indexed if taking it from # topology_reconstruction.get_all_current_points(). curr_points = topology_reconstruction.get_active_current_points() curr_points_including_inactive = ( topology_reconstruction.get_all_current_points() ) # Collect latitudes and longitudes of currently ACTIVE points in the ocean basin curr_lat_lon_points = [point.to_lat_lon() for point in curr_points] if curr_lat_lon_points: # Get the number of active points at this timestep. num_current_points = len(curr_points) # ndarray to fill with active point lats, lons and zvalues # FOR NOW, the number of gridding input columns is 6: # 0 = longitude # 1 = latitude # 2 = seafloor age # 3 = birth latitude snapshot # 4 = point id # 5 for the default gridding columns above, plus additional zvalues added next total_number_of_columns = 5 + len(self.zval_names) gridding_input_data = np.empty( [num_current_points, total_number_of_columns] ) # Lons and lats are first and second columns of the ndarray respectively gridding_input_data[:, 1], gridding_input_data[:, 0] = zip( *curr_lat_lon_points ) # NOTE: We need a single index to access data from curr_points_including_inactive AND allocate # this data to an ndarray with a number of rows equal to num_current_points. This index will # append +1 after each loop through curr_points_including_inactive. i = 0 # Get indices and points of all points at `time`, both active and inactive (which are NoneType points that # have undergone continental collision or subduction at `time`). for point_index, current_point in enumerate( curr_points_including_inactive ): # Look at all active points (these have not collided with a continent or trench) if current_point is not None: # Seafloor age gridding_input_data[i, 2] = ( appearance_time[point_index] - topology_reconstruction.get_current_time() ) # Birth latitude (snapshot) gridding_input_data[i, 3] = birth_lat[point_index] # Point ID (snapshot) gridding_input_data[i, 4] = point_id[ point_index ] # The ID of a corresponding point from the original list of all MOR-resolved points # GENERAL Z-VALUE ALLOCATION # Z values are 1st index onwards; 0th belongs to the point feature ID (thus +1) for j in range(len(self.zval_names)): # Adjusted index - and we have to add j to 5 to account for lat, lon, age, birth lat and point ID, adjusted_index = 5 + j # Spreading rate would be first # Access current zval from the master list of all zvalues for all points that ever existed in time_array gridding_input_data[i, adjusted_index] = zvalues[ point_index, j ] # Go to the next active point i += 1 gridding_input_dictionary = {} for i in list(range(total_number_of_columns)): gridding_input_dictionary[self.total_column_headers[i]] = [ list(j) for j in zip(*gridding_input_data) ][i] data_to_store = [ gridding_input_dictionary[i] for i in gridding_input_dictionary ] gridding_input_basename = "gridding_input_{:0.1f}Ma".format( topology_reconstruction.get_current_time() ) if self.file_collection is not None: gridding_input_basename = "{}_{}".format( self.file_collection, gridding_input_basename, ) gridding_input_filename = os.path.join( self.save_directory, gridding_input_basename ) # save debug file if ( "GPLATELY_DEBUG" in os.environ and os.environ["GPLATELY_DEBUG"].lower() == "true" ): save_age_grid_sample_points_to_gpml( gridding_input_dictionary["CURRENT_LONGITUDES"], gridding_input_dictionary["CURRENT_LATITUDES"], gridding_input_dictionary["SEAFLOOR_AGE"], topology_reconstruction.get_current_time(), self.save_directory, ) np.savez_compressed(gridding_input_filename, *data_to_store) if not topology_reconstruction.reconstruct_to_next_time(): break logger.info( f"Reconstruction done for {topology_reconstruction.get_current_time()}!" ) # return reconstruction_data def lat_lon_z_to_netCDF( self, zval_name, time_arr=None, unmasked=False, nprocs=1, ): """Produce a netCDF4 grid of a z-value identified by its `zval_name` for a given time range in `time_arr`. Seafloor age can be gridded by passing `zval_name` as `SEAFLOOR_AGE`, and spreading rate can be gridded with `SPREADING_RATE`. Saves all grids to compressed netCDF format in the attributed directory. Grids can be read into ndarray format using `gplately.grids.read_netcdf_grid()`. Parameters ---------- zval_name : str A string identifiers for a column in the ReconstructByTopologies gridding input files. time_arr : list of float, default None A time range to turn lons, lats and z-values into netCDF4 grids. If not provided, `time_arr` defaults to the full `time_array` provided to `SeafloorGrids`. unmasked : bool, default False Save unmasked grids, in addition to masked versions. nprocs : int, defaullt 1 Number of processes to use for certain operations (requires joblib). Passed to `joblib.Parallel`, so -1 means all available processes. """ parallel = None nprocs = int(nprocs) if nprocs != 1: try: from joblib import Parallel parallel = Parallel(nprocs) except ImportError: warnings.warn( "Could not import joblib; falling back to serial execution" ) # User can put any time array within SeafloorGrid bounds, but if none # is provided, it defaults to the attributed time array if time_arr is None: time_arr = self.time_array if parallel is None: for time in time_arr: _lat_lon_z_to_netCDF_time( time=time, zval_name=zval_name, file_collection=self.file_collection, save_directory=self.save_directory, total_column_headers=self.total_column_headers, extent=self.extent, resX=self.spacingX, resY=self.spacingY, unmasked=unmasked, continent_mask_filename=self.continent_mask_filename, ) else: from joblib import delayed parallel( delayed(_lat_lon_z_to_netCDF_time)( time=time, zval_name=zval_name, file_collection=self.file_collection, save_directory=self.save_directory, total_column_headers=self.total_column_headers, extent=self.extent, resX=self.spacingX, resY=self.spacingY, unmasked=unmasked, continent_mask_filename=self.continent_mask_filename, ) for time in time_arr )Instance variables
var PlotTopologiesTime-
Expand source code
@property def PlotTopologiesTime(self): return self._PlotTopologies_object.time var max_time-
The reconstruction time.
Expand source code
@property def max_time(self): """The reconstruction time.""" return self._max_time
Methods
def build_all_MOR_seedpoints(self)-
Resolve mid-ocean ridges for all times between
min_timeandmax_time, divide them into points that make up their shared sub-segments. Rotate these points to the left and right of the ridge using their stage rotation so that they spread from the ridge.Z-value allocation to each point is done here. In future, a function (like the spreading rate function) to calculate general z-data will be an input parameter.
Notes
If MOR seed point building is interrupted, progress is safeguarded as long as
resume_from_checkpointsis set toTrue.This assumes that points spread from ridges symmetrically, with the exception of large ridge jumps at successive timesteps. Therefore, z-values allocated to ridge-emerging points will appear symmetrical until changes in spreading ridge geometries create asymmetries.
In future, this will have a checkpoint save feature so that execution (which occurs during preparation for ReconstructByTopologies and can take several hours) can be safeguarded against run interruptions.
Parameters
time:float- The time at which to resolve ridge points and stage-rotate them off the ridge.
Returns
mor_point_features:FeatureCollection- All ridge seed points that have emerged from all ridge topologies at
time. These points have spread by being slightly rotated away from ridge locations attime.
References
get_mid_ocean_ridge_seedpoints() has been adapted from https://github.com/siwill22/agegrid-0.1/blob/master/automatic_age_grid_seeding.py#L117.
Expand source code
def build_all_MOR_seedpoints(self): """Resolve mid-ocean ridges for all times between `min_time` and `max_time`, divide them into points that make up their shared sub-segments. Rotate these points to the left and right of the ridge using their stage rotation so that they spread from the ridge. Z-value allocation to each point is done here. In future, a function (like the spreading rate function) to calculate general z-data will be an input parameter. Notes ----- If MOR seed point building is interrupted, progress is safeguarded as long as `resume_from_checkpoints` is set to `True`. This assumes that points spread from ridges symmetrically, with the exception of large ridge jumps at successive timesteps. Therefore, z-values allocated to ridge-emerging points will appear symmetrical until changes in spreading ridge geometries create asymmetries. In future, this will have a checkpoint save feature so that execution (which occurs during preparation for ReconstructByTopologies and can take several hours) can be safeguarded against run interruptions. Parameters ---------- time : float The time at which to resolve ridge points and stage-rotate them off the ridge. Returns ------- mor_point_features : FeatureCollection All ridge seed points that have emerged from all ridge topologies at `time`. These points have spread by being slightly rotated away from ridge locations at `time`. References ---------- get_mid_ocean_ridge_seedpoints() has been adapted from https://github.com/siwill22/agegrid-0.1/blob/master/automatic_age_grid_seeding.py#L117. """ # If we mustn't overwrite existing files in the `save_directory`, check the status of MOR seeding # to know where to start/continue seeding if self.resume_from_checkpoints: # Check the last MOR seedpoint gpmlz file that was built checkpointed_MOR_seedpoints = [ s.split("/")[-1] for s in glob.glob(self.save_directory + "/" + "*MOR_plus_one_points*") ] try: # -2 as an index accesses the age (float type), safeguards against identifying numbers in the SeafloorGrid.file_collection string last_seed_time = np.sort( [ float(re.findall(r"\d+", s)[-2]) for s in checkpointed_MOR_seedpoints ] )[0] # If none were built yet except: last_seed_time = "nil" # If MOR seeding has not started, start it from the top if last_seed_time == "nil": time_array = self.time_array # If the last seed time it could identify is outside the time bounds of the current instance of SeafloorGrid, start # from the top (this may happen if we use the same save directory for grids for a new set of times) elif last_seed_time not in self.time_array: time_array = self.time_array # If seeding was done to the min_time, we are finished elif last_seed_time == self.min_time: return # If seeding to `min_time` has been interrupted, resume it at last_masked_time. else: time_array = np.arange( last_seed_time, self.min_time - 0.1, -self.ridge_time_step ) # If we must overwrite all files in `save_directory`, start from `max_time`. else: time_array = self.time_array # Build all continental masks and spreading ridge points (with z values) self._get_mid_ocean_ridge_seedpoints(time_array) return def build_all_continental_masks(self)-
Create a continental mask to define the ocean basin for all times between
min_timeandmax_time. as well as to use as continental collision boundaries inReconstructByTopologies.Notes
Continental masking progress is safeguarded if ever masking is interrupted, provided that
resume_from_checkpointsis set toTrue.If
ReconstructByTopologiesidentifies a continental collision between oceanic points and the boundaries of this continental mask attime, those points are deleted attime.The continental mask is also saved to "/continent_mask_{}Ma.nc" as a compressed netCDF4 file if a
save_directoryis passed. Otherwise, the final grid is returned as a NumPy ndarray object.Returns
all_continental_masks:listofndarray- A masked grid per timestep in
time_arraywith 1=continental point, and 0=ocean point, for all points on the full global icosahedral mesh.
Expand source code
def build_all_continental_masks(self): """Create a continental mask to define the ocean basin for all times between `min_time` and `max_time`. as well as to use as continental collision boundaries in `ReconstructByTopologies`. Notes ----- Continental masking progress is safeguarded if ever masking is interrupted, provided that `resume_from_checkpoints` is set to `True`. If `ReconstructByTopologies` identifies a continental collision between oceanic points and the boundaries of this continental mask at `time`, those points are deleted at `time`. The continental mask is also saved to "/continent_mask_{}Ma.nc" as a compressed netCDF4 file if a `save_directory` is passed. Otherwise, the final grid is returned as a NumPy ndarray object. Returns ------- all_continental_masks : list of ndarray A masked grid per timestep in `time_array` with 1=continental point, and 0=ocean point, for all points on the full global icosahedral mesh. """ # If we mustn't overwrite existing files in the `save_directory`, check the status # of continental masking to know where to start/continue masking if self.resume_from_checkpoints: # Check the last continental mask that could be built checkpointed_continental_masks = [ s.split("/")[-1] for s in glob.glob(self.save_directory + "/" + "*continent_mask*") ] try: # -2 as an index accesses the age (float type), safeguards against identifying numbers in the SeafloorGrid.file_collection string last_masked_time = np.sort( [ float(re.findall(r"\d+", s)[-2]) for s in checkpointed_continental_masks ] )[0] # If none were built yet except: last_masked_time = "nil" # If masking has not started, start it from the top if last_masked_time == "nil": time_array = self.time_array # If the last seed time it could identify is outside the time bounds of the current instance of SeafloorGrid, start # from the top (this may happen if we use the same save directory for grids for a new set of times) elif last_masked_time not in self.time_array: time_array = self.time_array # If masking was done to the min_time, we are finished elif last_masked_time == self.min_time: return # If masking to `min_time` has been interrupted, resume it at last_masked_time. else: time_array = np.arange( last_masked_time, self.min_time - 0.1, -self.ridge_time_step ) # If we must overwrite all files in `save_directory`, start from `max_time`. else: time_array = self.time_array # Build all continental masks and spreading ridge points (with z values) self._create_continental_mask(time_array) return def create_initial_ocean_seed_points(self)-
Create the initial ocean basin seed point domain (at
max_timeonly) using Stripy's icosahedral triangulation with the specifiedself.refinement_levels.The ocean mesh starts off as a global-spanning Stripy icosahedral mesh.
create_initial_ocean_seed_pointspasses the automatically-resolved-to-current-time continental polygons from thePlotTopologies_object'scontinentsattribute (which can be from a COB terrane file or a continental polygon file) into Plate Tectonic Tools' point-in-polygon routine. It identifies ocean basin points that lie: * outside the polygons (for the ocean basin point domain) * inside the polygons (for the continental mask)Points from the mesh outside the continental polygons make up the ocean basin seed point mesh. The masked mesh is outputted as a compressed GPML (GPMLZ) file with the filename: "ocean_basin_seed_points_{}Ma.gpmlz" if a
save_directoryis passed. Otherwise, the mesh is returned as a pyGPlates FeatureCollection object.Notes
This point mesh represents ocean basin seafloor that was produced before
SeafloorGrid.max_time, and thus has unknown properties like valid time and spreading rate. As time passes, the plate reconstruction model sees points emerging from MORs. These new points spread to occupy the ocean basins, moving the initial filler points closer to subduction zones and continental polygons with which they can collide. If a collision is detected byPlateReconstructionsReconstructByTopologiesobject, these points are deleted.Ideally, if a reconstruction tree spans a large time range, all initial mesh points would collide with a continent or be subducted, leaving behind a mesh of well-defined MOR-emerged ocean basin points that data can be attributed to. However, some of these initial points situated close to contiental boundaries are retained through time - these form point artefacts with anomalously high ages. Even deep-time plate models (e.g. 1 Ga) will have these artefacts - removing them would require more detail to be added to the reconstruction model.
Returns
ocean_basin_point_mesh:pygplates.FeatureCollectionofpygplates.MultiPointOnSphere- A feature collection of point objects on the ocean basin.
Expand source code
def create_initial_ocean_seed_points(self): """Create the initial ocean basin seed point domain (at `max_time` only) using Stripy's icosahedral triangulation with the specified `self.refinement_levels`. The ocean mesh starts off as a global-spanning Stripy icosahedral mesh. `create_initial_ocean_seed_points` passes the automatically-resolved-to-current-time continental polygons from the `PlotTopologies_object`'s `continents` attribute (which can be from a COB terrane file or a continental polygon file) into Plate Tectonic Tools' point-in-polygon routine. It identifies ocean basin points that lie: * outside the polygons (for the ocean basin point domain) * inside the polygons (for the continental mask) Points from the mesh outside the continental polygons make up the ocean basin seed point mesh. The masked mesh is outputted as a compressed GPML (GPMLZ) file with the filename: "ocean_basin_seed_points_{}Ma.gpmlz" if a `save_directory` is passed. Otherwise, the mesh is returned as a pyGPlates FeatureCollection object. Notes ----- This point mesh represents ocean basin seafloor that was produced before `SeafloorGrid.max_time`, and thus has unknown properties like valid time and spreading rate. As time passes, the plate reconstruction model sees points emerging from MORs. These new points spread to occupy the ocean basins, moving the initial filler points closer to subduction zones and continental polygons with which they can collide. If a collision is detected by `PlateReconstruction`s `ReconstructByTopologies` object, these points are deleted. Ideally, if a reconstruction tree spans a large time range, **all** initial mesh points would collide with a continent or be subducted, leaving behind a mesh of well-defined MOR-emerged ocean basin points that data can be attributed to. However, some of these initial points situated close to contiental boundaries are retained through time - these form point artefacts with anomalously high ages. Even deep-time plate models (e.g. 1 Ga) will have these artefacts - removing them would require more detail to be added to the reconstruction model. Returns ------- ocean_basin_point_mesh : pygplates.FeatureCollection of pygplates.MultiPointOnSphere A feature collection of point objects on the ocean basin. """ if self.continent_mask_filename is None: # Ensure COB terranes at max time have reconstruction IDs and valid times COB_polygons = ensure_polygon_geometry( self._PlotTopologies_object.continents, self.rotation_model, self._max_time, ) # zval is a binary array encoding whether a point # coordinate is within a COB terrane polygon or not. # Use the icosahedral mesh MultiPointOnSphere attribute _, ocean_basin_point_mesh, zvals = point_in_polygon_routine( self.icosahedral_multi_point, COB_polygons ) # Plates to partition with plate_partitioner = pygplates.PlatePartitioner( COB_polygons, self.rotation_model, ) # Plate partition the ocean basin points meshnode_feature = pygplates.Feature( pygplates.FeatureType.create_from_qualified_string("gpml:MeshNode") ) meshnode_feature.set_geometry( ocean_basin_point_mesh # multi_point ) ocean_basin_meshnode = pygplates.FeatureCollection(meshnode_feature) paleogeography = plate_partitioner.partition_features( ocean_basin_meshnode, partition_return=pygplates.PartitionReturn.separate_partitioned_and_unpartitioned, properties_to_copy=[pygplates.PropertyName.gpml_shapefile_attributes], ) ocean_points = paleogeography[1] # Separate those inside polygons continent_points = paleogeography[0] # Separate those outside polygons # If a set of continent masks was passed, we can use max_time's continental # mask to build the initial profile of seafloor age. else: max_time_cont_mask = grids.Raster( self.continent_mask_filename.format(self._max_time) ) # If the input grid is at 0.5 degree uniform spacing, then the input # grid is 7x more populated than a 6-level stripy icosahedral mesh and # using this resolution for the initial ocean mesh will dramatically slow down # reconstruction by topologies. # Scale down the resolution based on the input mask resolution # (percentage was found in __init__.) max_time_cont_mask.resize( int(max_time_cont_mask.shape[0] * self.percentage), int(max_time_cont_mask.shape[1] * self.percentage), inplace=True, ) lat = np.linspace(-90, 90, max_time_cont_mask.shape[0]) lon = np.linspace(-180, 180, max_time_cont_mask.shape[1]) llon, llat = np.meshgrid(lon, lat) mask_inds = np.where(max_time_cont_mask.data.flatten() == 0) mask_vals = max_time_cont_mask.data.flatten() mask_lon = llon.flatten()[mask_inds] mask_lat = llat.flatten()[mask_inds] ocean_pt_feature = pygplates.Feature() ocean_pt_feature.set_geometry( pygplates.MultiPointOnSphere(zip(mask_lat, mask_lon)) ) ocean_points = [ocean_pt_feature] # Now that we have ocean points... # Determine age of ocean basin points using their proximity to MOR features # and an assumed globally-uniform ocean basin mean spreading rate. # We need resolved topologies at the `max_time` to pass into the proximity # function resolved_topologies = [] shared_boundary_sections = [] pygplates.resolve_topologies( self.topology_features, self.rotation_model, resolved_topologies, self._max_time, shared_boundary_sections, ) pX, pY, pZ = tools.find_distance_to_nearest_ridge( resolved_topologies, shared_boundary_sections, ocean_points, ) # Divide spreading rate by 2 to use half the mean spreading rate pAge = np.array(pZ) / (self.initial_ocean_mean_spreading_rate / 2.0) initial_ocean_point_features = [] initial_ocean_multipoints = [] for point in zip(pX, pY, pAge): point_feature = pygplates.Feature() point_feature.set_geometry(pygplates.PointOnSphere(point[1], point[0])) # Add 'time' to the age at the time of computation, to get the valid time in Ma point_feature.set_valid_time(point[2] + self._max_time, -1) # For now: custom zvals are added as shapefile attributes - will attempt pandas data frames # point_feature = set_shapefile_attribute(point_feature, self.initial_ocean_mean_spreading_rate, "SPREADING_RATE") # Seems like static data initial_ocean_point_features.append(point_feature) initial_ocean_multipoints.append(point_feature.get_geometry()) # print(initial_ocean_point_features) multi_point_feature = pygplates.MultiPointOnSphere(initial_ocean_multipoints) basename = "ocean_basin_seed_points_{}_RLs_{}Ma.gpmlz".format( self.refinement_levels, self._max_time, ) if self.file_collection is not None: basename = "{}_{}".format(self.file_collection, basename) output_filename = os.path.join(self.save_directory, basename) initial_ocean_feature_collection = pygplates.FeatureCollection( initial_ocean_point_features ) initial_ocean_feature_collection.write(output_filename) # Collect all point feature data into a pandas dataframe self._collect_point_data_in_dataframe( initial_ocean_feature_collection, np.array( [self.initial_ocean_mean_spreading_rate] * len(pX) ), # for now, spreading rate is one zvalue for initial ocean points. will other zvalues need to have a generalised workflow? self._max_time, ) return ( pygplates.FeatureCollection(initial_ocean_point_features), multi_point_feature, ) def lat_lon_z_to_netCDF(self, zval_name, time_arr=None, unmasked=False, nprocs=1)-
Produce a netCDF4 grid of a z-value identified by its
zval_namefor a given time range intime_arr.Seafloor age can be gridded by passing
zval_nameasSEAFLOOR_AGE, and spreading rate can be gridded withSPREADING_RATE.Saves all grids to compressed netCDF format in the attributed directory. Grids can be read into ndarray format using
read_netcdf_grid().Parameters
zval_name:str- A string identifiers for a column in the ReconstructByTopologies gridding input files.
time_arr:listoffloat, defaultNone- A time range to turn lons, lats and z-values into netCDF4 grids. If not provided,
time_arrdefaults to the fulltime_arrayprovided toSeafloorGrids. unmasked:bool, defaultFalse- Save unmasked grids, in addition to masked versions.
nprocs:int, defaullt 1- Number of processes to use for certain operations (requires joblib).
Passed to
joblib.Parallel, so -1 means all available processes.
Expand source code
def lat_lon_z_to_netCDF( self, zval_name, time_arr=None, unmasked=False, nprocs=1, ): """Produce a netCDF4 grid of a z-value identified by its `zval_name` for a given time range in `time_arr`. Seafloor age can be gridded by passing `zval_name` as `SEAFLOOR_AGE`, and spreading rate can be gridded with `SPREADING_RATE`. Saves all grids to compressed netCDF format in the attributed directory. Grids can be read into ndarray format using `gplately.grids.read_netcdf_grid()`. Parameters ---------- zval_name : str A string identifiers for a column in the ReconstructByTopologies gridding input files. time_arr : list of float, default None A time range to turn lons, lats and z-values into netCDF4 grids. If not provided, `time_arr` defaults to the full `time_array` provided to `SeafloorGrids`. unmasked : bool, default False Save unmasked grids, in addition to masked versions. nprocs : int, defaullt 1 Number of processes to use for certain operations (requires joblib). Passed to `joblib.Parallel`, so -1 means all available processes. """ parallel = None nprocs = int(nprocs) if nprocs != 1: try: from joblib import Parallel parallel = Parallel(nprocs) except ImportError: warnings.warn( "Could not import joblib; falling back to serial execution" ) # User can put any time array within SeafloorGrid bounds, but if none # is provided, it defaults to the attributed time array if time_arr is None: time_arr = self.time_array if parallel is None: for time in time_arr: _lat_lon_z_to_netCDF_time( time=time, zval_name=zval_name, file_collection=self.file_collection, save_directory=self.save_directory, total_column_headers=self.total_column_headers, extent=self.extent, resX=self.spacingX, resY=self.spacingY, unmasked=unmasked, continent_mask_filename=self.continent_mask_filename, ) else: from joblib import delayed parallel( delayed(_lat_lon_z_to_netCDF_time)( time=time, zval_name=zval_name, file_collection=self.file_collection, save_directory=self.save_directory, total_column_headers=self.total_column_headers, extent=self.extent, resX=self.spacingX, resY=self.spacingY, unmasked=unmasked, continent_mask_filename=self.continent_mask_filename, ) for time in time_arr ) def prepare_for_reconstruction_by_topologies(self)-
Prepare three main auxiliary files for seafloor data gridding: * Initial ocean seed points (at
max_time) * Continental masks (frommax_timetomin_time) * MOR points (frommax_timetomin_time)Returns lists of all attributes for the initial ocean point mesh and all ridge points for all times in the reconstruction time array.
Expand source code
def prepare_for_reconstruction_by_topologies(self): """Prepare three main auxiliary files for seafloor data gridding: * Initial ocean seed points (at `max_time`) * Continental masks (from `max_time` to `min_time`) * MOR points (from `max_time` to `min_time`) Returns lists of all attributes for the initial ocean point mesh and all ridge points for all times in the reconstruction time array. """ # INITIAL OCEAN SEED POINT MESH ---------------------------------------------------- ( initial_ocean_seed_points, initial_ocean_seed_points_mp, ) = self.create_initial_ocean_seed_points() logger.info("Finished building initial_ocean_seed_points!") # MOR SEED POINTS AND CONTINENTAL MASKS -------------------------------------------- # The start time for seeding is controlled by whether the overwrite_existing_gridding_inputs # parameter is set to `True` (in which case the start time is `max_time`). If it is `False` # and; # - a run of seeding and continental masking was interrupted, and ridge points were # checkpointed at n Ma, seeding resumes at n-1 Ma until `min_time` or another interruption # occurs; # - seeding was completed but the subsequent gridding input creation was interrupted, # seeding is assumed completed and skipped. The workflow automatically proceeds to re-gridding. if self.continent_mask_filename is None: self.build_all_continental_masks() else: logger.info( "Continent masks passed to SeafloorGrid - skipping continental mask generation!" ) self.build_all_MOR_seedpoints() # ALL-TIME POINTS ----------------------------------------------------- # Extract all feature attributes for all reconstruction times into lists active_points = [] appearance_time = [] birth_lat = [] # latitude_of_crust_formation prev_lat = [] prev_lon = [] # Extract point feature attributes from MOR seed points all_mor_features = [] zvalues = np.empty((0, len(self.zval_names))) for time in self.time_array: # If we're at the maximum time, start preparing points from the initial ocean mesh # as well as their z values if time == self._max_time: for feature in initial_ocean_seed_points: active_points.append(feature.get_geometry()) appearance_time.append(feature.get_valid_time()[0]) birth_lat.append(feature.get_geometry().to_lat_lon_list()[0][0]) prev_lat.append(feature.get_geometry().to_lat_lon_list()[0][0]) prev_lon.append(feature.get_geometry().to_lat_lon_list()[0][1]) curr_zvalues = self._extract_zvalues_from_npz_to_ndarray( initial_ocean_seed_points, time ) zvalues = np.concatenate((zvalues, curr_zvalues), axis=0) # Otherwise, we'd be preparing MOR points and their z values else: # GPMLZ file of MOR seedpoints basename = "MOR_plus_one_points_{:0.2f}.gpmlz".format(time) if self.file_collection is not None: basename = "{}_{}".format(self.file_collection, basename) filename = os.path.join(self.save_directory, basename) features = pygplates.FeatureCollection(filename) for feature in features: if feature.get_valid_time()[0] < self.time_array[0]: active_points.append(feature.get_geometry()) appearance_time.append(feature.get_valid_time()[0]) birth_lat.append(feature.get_geometry().to_lat_lon_list()[0][0]) prev_lat.append(feature.get_geometry().to_lat_lon_list()[0][0]) prev_lon.append(feature.get_geometry().to_lat_lon_list()[0][1]) # COLLECT NDARRAY OF ALL ZVALUES IN THIS TIMESTEP ------------------ curr_zvalues = self._extract_zvalues_from_npz_to_ndarray(features, time) zvalues = np.concatenate((zvalues, curr_zvalues), axis=0) return active_points, appearance_time, birth_lat, prev_lat, prev_lon, zvalues def reconstruct_by_topologies(self)-
Obtain all active ocean seed points at
time- these are points that have not been consumed at subduction zones or have not collided with continental polygons.All active points' latitudes, longitues, seafloor ages, spreading rates and all other general z-values are saved to a gridding input file (.npz).
Expand source code
def reconstruct_by_topologies(self): """Obtain all active ocean seed points at `time` - these are points that have not been consumed at subduction zones or have not collided with continental polygons. All active points' latitudes, longitues, seafloor ages, spreading rates and all other general z-values are saved to a gridding input file (.npz). """ logger.info("Preparing all initial files...") # Obtain all info from the ocean seed points and all MOR points through time, store in # arrays ( active_points, appearance_time, birth_lat, prev_lat, prev_lon, zvalues, ) = self.prepare_for_reconstruction_by_topologies() #### Begin reconstruction by topology process: # Indices for all points (`active_points`) that have existed from `max_time` to `min_time`. point_id = range(len(active_points)) # Specify the default collision detection region as subduction zones default_collision = reconstruction._DefaultCollision( feature_specific_collision_parameters=[ ( pygplates.FeatureType.gpml_subduction_zone, self.subduction_collision_parameters, ) ] ) # In addition to the default subduction detection, also detect continental collisions # Use the input continent mask if it is provided. if self.continent_mask_filename is not None: collision_spec = reconstruction._ContinentCollision( # This filename string should not have a time formatted into it - this is # taken care of later. self.continent_mask_filename, default_collision, verbose=False, ) else: # If a continent mask is not provided, use the ones made. mask_basename = r"continent_mask_{}Ma.nc" if self.file_collection is not None: mask_basename = str(self.file_collection) + "_" + mask_basename mask_template = os.path.join(self.save_directory, mask_basename) collision_spec = reconstruction._ContinentCollision( mask_template, default_collision, verbose=False, ) # Call the reconstruct by topologies object topology_reconstruction = reconstruction._ReconstructByTopologies( self.rotation_model, self.topology_features, self._max_time, self.min_time, self.ridge_time_step, active_points, point_begin_times=appearance_time, detect_collisions=collision_spec, ) # Initialise the reconstruction. topology_reconstruction.begin_reconstruction() # Loop over the reconstruction times until the end of the reconstruction time span, or until # all points have entered their valid time range *and* either exited their time range or # have been deactivated (subducted forward in time or consumed by MOR backward in time). reconstruction_data = [] while True: logger.info( f"Reconstruct by topologies: working on time {topology_reconstruction.get_current_time():0.2f} Ma" ) # NOTE: # topology_reconstruction.get_active_current_points() and topology_reconstruction.get_all_current_points() # are different. The former is a subset of the latter, and it represents all points at the timestep that # have not collided with a continental or subduction boundary. The remainders in the latter are inactive # (NoneType) points, which represent the collided points. # We need to access active point data from topology_reconstruction.get_all_current_points() because it has # the same length as the list of all initial ocean points and MOR seed points that have ever emerged from # spreading ridge topologies through `max_time` to `min_time`. Therefore, it protects the time and space # order in which all MOR points through time were seeded by pyGPlates. At any given timestep, not all these # points will be active, but their indices are retained. Thus, z value allocation, point latitudes and # longitudes of active points will be correctly indexed if taking it from # topology_reconstruction.get_all_current_points(). curr_points = topology_reconstruction.get_active_current_points() curr_points_including_inactive = ( topology_reconstruction.get_all_current_points() ) # Collect latitudes and longitudes of currently ACTIVE points in the ocean basin curr_lat_lon_points = [point.to_lat_lon() for point in curr_points] if curr_lat_lon_points: # Get the number of active points at this timestep. num_current_points = len(curr_points) # ndarray to fill with active point lats, lons and zvalues # FOR NOW, the number of gridding input columns is 6: # 0 = longitude # 1 = latitude # 2 = seafloor age # 3 = birth latitude snapshot # 4 = point id # 5 for the default gridding columns above, plus additional zvalues added next total_number_of_columns = 5 + len(self.zval_names) gridding_input_data = np.empty( [num_current_points, total_number_of_columns] ) # Lons and lats are first and second columns of the ndarray respectively gridding_input_data[:, 1], gridding_input_data[:, 0] = zip( *curr_lat_lon_points ) # NOTE: We need a single index to access data from curr_points_including_inactive AND allocate # this data to an ndarray with a number of rows equal to num_current_points. This index will # append +1 after each loop through curr_points_including_inactive. i = 0 # Get indices and points of all points at `time`, both active and inactive (which are NoneType points that # have undergone continental collision or subduction at `time`). for point_index, current_point in enumerate( curr_points_including_inactive ): # Look at all active points (these have not collided with a continent or trench) if current_point is not None: # Seafloor age gridding_input_data[i, 2] = ( appearance_time[point_index] - topology_reconstruction.get_current_time() ) # Birth latitude (snapshot) gridding_input_data[i, 3] = birth_lat[point_index] # Point ID (snapshot) gridding_input_data[i, 4] = point_id[ point_index ] # The ID of a corresponding point from the original list of all MOR-resolved points # GENERAL Z-VALUE ALLOCATION # Z values are 1st index onwards; 0th belongs to the point feature ID (thus +1) for j in range(len(self.zval_names)): # Adjusted index - and we have to add j to 5 to account for lat, lon, age, birth lat and point ID, adjusted_index = 5 + j # Spreading rate would be first # Access current zval from the master list of all zvalues for all points that ever existed in time_array gridding_input_data[i, adjusted_index] = zvalues[ point_index, j ] # Go to the next active point i += 1 gridding_input_dictionary = {} for i in list(range(total_number_of_columns)): gridding_input_dictionary[self.total_column_headers[i]] = [ list(j) for j in zip(*gridding_input_data) ][i] data_to_store = [ gridding_input_dictionary[i] for i in gridding_input_dictionary ] gridding_input_basename = "gridding_input_{:0.1f}Ma".format( topology_reconstruction.get_current_time() ) if self.file_collection is not None: gridding_input_basename = "{}_{}".format( self.file_collection, gridding_input_basename, ) gridding_input_filename = os.path.join( self.save_directory, gridding_input_basename ) # save debug file if ( "GPLATELY_DEBUG" in os.environ and os.environ["GPLATELY_DEBUG"].lower() == "true" ): save_age_grid_sample_points_to_gpml( gridding_input_dictionary["CURRENT_LONGITUDES"], gridding_input_dictionary["CURRENT_LATITUDES"], gridding_input_dictionary["SEAFLOOR_AGE"], topology_reconstruction.get_current_time(), self.save_directory, ) np.savez_compressed(gridding_input_filename, *data_to_store) if not topology_reconstruction.reconstruct_to_next_time(): break logger.info( f"Reconstruction done for {topology_reconstruction.get_current_time()}!" ) # return reconstruction_data def update_time(self, max_time)-
Expand source code
def update_time(self, max_time): self._max_time = float(max_time) self._PlotTopologies_object.time = float(max_time)