Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Properly implement non-zero anchor plates #131

Merged
merged 3 commits into from
Nov 9, 2023
Merged

Properly implement non-zero anchor plates #131

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e7a9ca1
add Band1 as a recognised z variable
Nov 9, 2023
6f8669d
fix anchor plate IDs for subduction zone, MOR tessellation functions
Nov 9, 2023
6e33a7b
implement in reconstruct and reconstruct_to_birth_age in Points object
Nov 9, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion gplately/grids.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -180,7 +180,7 @@ def find_label(keys, labels):
# possible permutations of lon/lat/z
label_lon = ['lon', 'lons', 'longitude', 'x', 'east', 'easting', 'eastings']
label_lat = ['lat', 'lats', 'latitude', 'y', 'north', 'northing', 'northings']
label_z = ['z', 'data', 'values']
label_z = ['z', 'data', 'values', 'Band1']

# add capitalise and upper case permutations
label_lon = label_lon + [label.capitalize() for label in label_lon] + [label.upper() for label in label_lon]
Expand Down
58 changes: 41 additions & 17 deletions gplately/reconstruction.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,21 +45,23 @@ def __init__(
rotation_model,
topology_features=None,
static_polygons=None,
default_anchor_plate_id=0,
anchor_plate_id=0,
):
if hasattr(rotation_model, "reconstruction_identifier"):
self.name = rotation_model.reconstruction_identifier
else:
self.name = None

self.anchor_plate_id = int(anchor_plate_id)
self.rotation_model = _RotationModel(
rotation_model, default_anchor_plate_id=default_anchor_plate_id
rotation_model, default_anchor_plate_id=anchor_plate_id
)
self.topology_features = _load_FeatureCollection(topology_features)
self.static_polygons = _load_FeatureCollection(static_polygons)

def __getstate__(self):
filenames = {"rotation_model": self.rotation_model.filenames}
filenames = {"rotation_model": self.rotation_model.filenames,
"anchor_plate_id": self.anchor_plate_id}

if self.topology_features:
filenames["topology_features"] = self.topology_features.filenames
Expand All @@ -78,7 +80,7 @@ def __getstate__(self):

def __setstate__(self, state):
# reinstate unpicklable items
self.rotation_model = _RotationModel(state["rotation_model"])
self.rotation_model = _RotationModel(state["rotation_model"], default_anchor_plate_id=state["anchor_plate_id"])

self.topology_features = None
self.static_polygons = None
Expand Down Expand Up @@ -172,6 +174,8 @@ def tessellate_subduction_zones(

The delta time interval used for velocity calculations is, by default, assumed to be 1Ma.
"""
anchor_plate_id = kwargs.pop("anchor_plate_id", self.anchor_plate_id)

if ignore_warnings:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
Expand All @@ -180,6 +184,7 @@ def tessellate_subduction_zones(
self.topology_features,
tessellation_threshold_radians,
float(time),
anchor_plate_id=anchor_plate_id,
**kwargs
)

Expand All @@ -189,6 +194,7 @@ def tessellate_subduction_zones(
self.topology_features,
tessellation_threshold_radians,
float(time),
anchor_plate_id=anchor_plate_id,
**kwargs
)

Expand Down Expand Up @@ -464,6 +470,8 @@ def tessellate_mid_ocean_ridges(
* spreading velocity magnitude (in cm/yr)
* length of arc segment (in degrees) that current point is on
"""
anchor_plate_id = kwargs.pop("anchor_plate_id", self.anchor_plate_id)

if ignore_warnings:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
Expand All @@ -474,6 +482,7 @@ def tessellate_mid_ocean_ridges(
float(time),
tessellation_threshold_radians,
spreading_feature_types,
anchor_plate_id=anchor_plate_id,
**kwargs
)

Expand All @@ -485,6 +494,7 @@ def tessellate_mid_ocean_ridges(
float(time),
tessellation_threshold_radians,
spreading_feature_types,
anchor_plate_id=anchor_plate_id,
**kwargs
)

Expand Down Expand Up @@ -599,7 +609,7 @@ def total_ridge_length(self, time, use_ptt=False, ignore_warnings=False):

return total_ridge_length_kms

def reconstruct(self, feature, to_time, from_time=0, anchor_plate_id=0, **kwargs):
def reconstruct(self, feature, to_time, from_time=0, anchor_plate_id=None, **kwargs):
"""Reconstructs regular geological features, motion paths or flowlines to a specific geological time.

Parameters
Expand All @@ -616,10 +626,10 @@ def reconstruct(self, feature, to_time, from_time=0, anchor_plate_id=0, **kwargs
The specific geological time to reconstruct from. By default, this is set to present day. Raises
`NotImplementedError` if `from_time` is not set to 0.0 Ma (present day).

anchor_plate_id : int, default=0
anchor_plate_id : int, default=None
Reconstruct features with respect to a certain anchor plate. By default, reconstructions are made
with respect to the absolute reference frame, like a stationary geological element (e.g. a mantle
plume). This frame is given the plate ID of 0.
with respect to the absolute reference frame (anchor_plate_id = 0), like a stationary object in the mantle,
unless otherwise specified.

**reconstruct_type : ReconstructType, default=ReconstructType.feature_geometry
The specific reconstruction type to generate based on input feature geometry type. Can be provided as
Expand Down Expand Up @@ -657,6 +667,10 @@ def reconstruct(self, feature, to_time, from_time=0, anchor_plate_id=0, **kwargs
from_time, to_time = float(from_time), float(to_time)

reconstructed_features = []

if not anchor_plate_id:
anchor_plate_id = self.anchor_plate_id

pygplates.reconstruct(
feature,
self.rotation_model,
Expand Down Expand Up @@ -780,7 +794,7 @@ def create_motion_path(
lats,
time_array,
plate_id=None,
anchor_plate_id=0,
anchor_plate_id=None,
return_rate_of_motion=False,
):
"""Create a path of points to mark the trajectory of a plate's motion
Expand Down Expand Up @@ -846,6 +860,9 @@ def create_motion_path(
query_plate_id = False
plate_ids = np.ones(len(lons), dtype=int) * plate_id

if not anchor_plate_id:
anchor_plate_id = self.anchor_plate_id

seed_points = zip(lats, lons)
for i, lat_lon in enumerate(seed_points):
seed_points_at_digitisation_time = pygplates.PointOnSphere(
Expand Down Expand Up @@ -1403,7 +1420,7 @@ def get_geodataframe(self):
"""
return self.get_geopandas_dataframe()

def reconstruct(self, time, anchor_plate_id=0, return_array=False, **kwargs):
def reconstruct(self, time, anchor_plate_id=None, return_array=False, **kwargs):
"""Reconstructs regular geological features, motion paths or flowlines to a specific geological time and extracts
the latitudinal and longitudinal points of these features.

Expand All @@ -1415,10 +1432,10 @@ def reconstruct(self, time, anchor_plate_id=0, return_array=False, **kwargs):
time : float
The specific geological time (Ma) to reconstruct features to.

anchor_plate_id : int, default=0
anchor_plate_id : int, default=None
Reconstruct features with respect to a certain anchor plate. By default, reconstructions are made
with respect to the absolute reference frame, like a stationary geological element (e.g. a mantle
plume). This frame is given the plate ID of 0.
with respect to the anchor_plate_ID specified in the `gplately.PlateReconstruction` object,
which is a default plate ID of 0 unless otherwise specified.

return_array : bool, default False
Return a `numpy.ndarray`, rather than a `Points` object.
Expand Down Expand Up @@ -1455,6 +1472,10 @@ def reconstruct(self, time, anchor_plate_id=0, return_array=False, **kwargs):
"""
from_time = self.time
to_time = time

if not anchor_plate_id:
anchor_plate_id = self.plate_reconstruction.anchor_plate_id

reconstructed_features = self.plate_reconstruction.reconstruct(
self.features, to_time, from_time, anchor_plate_id=anchor_plate_id, **kwargs
)
Expand All @@ -1473,7 +1494,7 @@ def reconstruct(self, time, anchor_plate_id=0, return_array=False, **kwargs):
gpts.add_attributes(**self.attributes.copy())
return gpts

def reconstruct_to_birth_age(self, ages, anchor_plate_id=0, **kwargs):
def reconstruct_to_birth_age(self, ages, anchor_plate_id=None, **kwargs):
"""Reconstructs point features supplied to the `Points` object from the supplied initial time (`self.time`)
to a range of times. The number of supplied times must equal the number of point features supplied to the Points object.

Expand All @@ -1482,10 +1503,10 @@ def reconstruct_to_birth_age(self, ages, anchor_plate_id=0, **kwargs):
ages : array
Geological times to reconstruct features to. Must have the same length as the `Points `object's `self.features` attribute
(which holds all point features represented on a unit length sphere in 3D Cartesian coordinates).
anchor_plate_id : int, default=0
anchor_plate_id : int, default=None
Reconstruct features with respect to a certain anchor plate. By default, reconstructions are made
with respect to the absolute reference frame, like a stationary geological element (e.g. a mantle
plume). This frame is given the plate ID of 0.
with respect to the anchor_plate_ID specified in the `gplately.PlateReconstruction` object,
which is a default plate ID of 0 unless otherwise specified.
**kwargs
Additional keyword arguments for the `gplately.PlateReconstruction.reconstruct` method.

Expand Down Expand Up @@ -1517,6 +1538,9 @@ def reconstruct_to_birth_age(self, ages, anchor_plate_id=0, **kwargs):

"""
from_time = self.time
if not anchor_plate_id:
anchor_plate_id = self.plate_reconstruction.anchor_plate_id

ages = np.array(ages)

if len(ages) != len(self.features):
Expand Down
Loading