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

Optionally change reference frames #3

Merged
merged 1 commit into from
Jul 25, 2024
Merged

Optionally change reference frames #3

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
88 changes: 76 additions & 12 deletions slabdip/predictor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,14 @@ def model(self):
def model(self, value):
self.add_plate_reconstruction(value)

def sample_age_grid(self, lons, lats, time):
def sample_age_grid(
self,
lons,
lats,
time,
from_rotation_reference_plate=None,
to_rotation_reference_plate=None
):
# age_grid = self.downloader.get_age_grid(reconstruction_time)
if self.agegrid_filename:
age_raster = gplately.Raster(data=self.agegrid_filename.format(time))
Expand All @@ -125,10 +132,32 @@ def sample_age_grid(self, lons, lats, time):
Set agegrid_filename or provide a supported reconstruction model")

age_raster.fill_NaNs(inplace=True) # fill in NaN values

# rotate grids if needed
if from_rotation_reference_plate is not None and to_rotation_reference_plate is not None:
spacingX = np.diff(age_raster.lons).mean()
spacingY = np.diff(age_raster.lats).mean()
mean_spacing = 0.5*(spacingX + spacingY)
age_raster = age_raster.rotate_reference_frames(
mean_spacing,
time,
from_rotation_features_or_model=self._model.rotation_model,
to_rotation_features_or_model=self._model.rotation_model,
from_rotation_reference_plate=from_rotation_reference_plate,
to_rotation_reference_plate=to_rotation_reference_plate,
)

age_interp = age_raster.interpolate(lons, lats) # interpolate to trenches
return age_interp

def sample_spreading_rate_grid(self, lons, lats, time):
def sample_spreading_rate_grid(
self,
lons,
lats,
time,
from_rotation_reference_plate=None,
to_rotation_reference_plate=None
):
# spreadrate_grid = self.downloader.get_spreading_rate_grid(reconstruction_time)
if self.spreadrate_filename:
spreadrate_raster = gplately.Raster(data=self.spreadrate_filename.format(time))
Expand All @@ -140,6 +169,21 @@ def sample_spreading_rate_grid(self, lons, lats, time):
Set spreadrate_filename or provide a supported reconstruction model")

spreadrate_raster.fill_NaNs(inplace=True)

# rotate grids if needed
if from_rotation_reference_plate is not None and to_rotation_reference_plate is not None:
spacingX = np.diff(spreadrate_raster.lons).mean()
spacingY = np.diff(spreadrate_raster.lats).mean()
mean_spacing = 0.5*(spacingX + spacingY)
spreadrate_raster = spreadrate_raster.rotate_reference_frames(
mean_spacing,
time,
from_rotation_features_or_model=self._model.rotation_model,
to_rotation_features_or_model=self._model.rotation_model,
from_rotation_reference_plate=from_rotation_reference_plate,
to_rotation_reference_plate=to_rotation_reference_plate,
)

spreadrate_interp = spreadrate_raster.interpolate(lons, lats)
return spreadrate_interp*1e-3

Expand Down Expand Up @@ -248,23 +292,29 @@ def calculate_trench_curvature(self, lons, lats, norm, length, smoothing=5, retu
else:
return subduction_radius

def tessellate_slab_dip(self, time, tessellation_threshold_radians=DEFAULT_TESSELLATION):
def tessellate_slab_dip(
self,
time,
tessellation_threshold_radians=DEFAULT_TESSELLATION,
from_rotation_reference_plate=None,
to_rotation_reference_plate=None
):
time = int(time)
self.tessellation_threshold_radians = DEFAULT_TESSELLATION
self.tessellation_threshold_radians = tessellation_threshold_radians

if self._model is None:
raise ValueError("Don't forget to set a GPlately plate model! `self.model = model`")


# update plate ID
if to_rotation_reference_plate is not None:
self._model.anchor_plate_id = to_rotation_reference_plate

subduction_data = self._model.tessellate_subduction_zones(
time,
tessellation_threshold_radians,
ignore_warnings=True,
output_subducting_absolute_velocity_components=True)

# mask "negative" subduction rates
subduction_convergence = np.fabs(subduction_data[:,2])*1e-2 * np.cos(np.deg2rad(subduction_data[:,3]))
subduction_data = subduction_data[subduction_convergence >= 0]

subduction_lon = subduction_data[:,0]
subduction_lat = subduction_data[:,1]
subduction_vel = subduction_data[:,2]*1e-2
Expand All @@ -277,12 +327,26 @@ def tessellate_slab_dip(self, time, tessellation_threshold_radians=DEFAULT_TESSE
subduction_migration = np.fabs(subduction_data[:,4])*1e-2 * np.cos(np.deg2rad(subduction_data[:,5]))
subduction_plate_vel = subduction_data[:,10]*1e-2

subduction_convergence = np.clip(subduction_convergence, 0, 1e99)

# sample AgeGrid
age_interp = self.sample_age_grid(subduction_lon, subduction_lat, time)
age_interp = self.sample_age_grid(
subduction_lon,
subduction_lat,
time,
from_rotation_reference_plate=from_rotation_reference_plate,
to_rotation_reference_plate=to_rotation_reference_plate)

# calculate the thickness of the downgoing plate
thickness = gplately.tools.plate_isotherm_depth(age_interp)

# sample spreadrate grid
spreadrate_interp = self.sample_spreading_rate_grid(subduction_lon, subduction_lat, time)
spreadrate_interp = self.sample_spreading_rate_grid(
subduction_lon,
subduction_lat,
time,
from_rotation_reference_plate=from_rotation_reference_plate,
to_rotation_reference_plate=to_rotation_reference_plate)

# get the ratio of convergence velocity to trench migration
vratio = (subduction_convergence + subduction_migration)/(subduction_convergence + 1e-22)
Expand Down