# Split an isochron into ridge and transform segments

This example splits the geometry of an isochron into ridge and transform segments based on each segment’s alignment with the isochron’s stage pole at its time of appearance.

Warning

This algorithm only works well under the following conditions:

ridge segments are perpendicular to their spreading directions,

isochron geometries are up-to-date with respect to the rotation model (ie, stage pole is in correct location relative to geometry),

there are valid rotations (in rotation model) for each isochron at its birth time plus one (ie, 1My

*prior*to isochron birth time),all isochrons have

`conjugate plate IDs`

.

Note that you can tweak the *isochron_segment_deviation_in_radians* parameter. A 45 degree split may not always be appropriate.

## Sample code

```
import pygplates
import math
# How much an isochron segment can deviate from the stage pole before it's considered a transform segment.
isochron_segment_deviation_in_radians = math.pi / 4 # An even 45 degrees split
# Load one or more rotation files into a rotation model.
rotation_model = pygplates.RotationModel('rotations.rot')
# Load the isochron features.
isochron_features = pygplates.FeatureCollection('isochrons.gpml')
# Iterate over all geometries in isochron features.
for isochron_feature in isochron_features:
begin_time, end_time = isochron_feature.get_valid_time()
plate_id = isochron_feature.get_reconstruction_plate_id()
conjugate_plate_id = isochron_feature.get_conjugate_plate_id()
# Calculate the stage rotation at the isochron birth time of the isochron's plate relative
# to its conjugate plate.
stage_rotation = rotation_model.get_rotation(
begin_time + 1, plate_id, begin_time, conjugate_plate_id, conjugate_plate_id)
stage_pole, stage_angle_radians = stage_rotation.get_euler_pole_and_angle()
# Present day geometries need to be rotated, relative to the conjugate plate, to the 'from' time
# of the above stage rotation (which is 'begin_time') so that they can then be rotated by
# the stage rotation. To avoid having to rotate the present day geometries into this stage pole
# reference frame we can instead apply the inverse rotation to the stage pole itself.
stage_pole_reference_frame = rotation_model.get_rotation(
begin_time, plate_id, 0, conjugate_plate_id, conjugate_plate_id)
stage_pole = stage_pole_reference_frame.get_inverse() * stage_pole
# A feature usually has a single geometry but it could have more - iterate over them all.
# Note that we are iterating over the un-rotated (or present day) geometries as noted above.
for isochron_geometry in isochron_feature.get_geometries():
# Group the current isochron geometry into ridge and transform segments.
ridge_segments = []
transform_segments = []
# Iterate over the segments of the current geometry.
# Note that we're assuming the geometry is a polyline (or polygon) - otherwise this will raise an error.
for segment in isochron_geometry.get_segments():
# Ignore zero length segments - they don't have a direction.
if segment.is_zero_length():
continue
# Get the point in the middle of the segment and its tangential direction.
segment_midpoint = segment.get_arc_point(0.5)
segment_direction_at_midpoint = segment.get_arc_direction(0.5)
# Get the direction from the segment midpoint to the stage pole.
# This is the tangential direction at the start of an arc from the segment
# midpoint to the stage pole (the zero parameter indicates the arc start point
# which is the segment midpoint).
segment_to_stage_pole_direction = pygplates.GreatCircleArc(
segment_midpoint, stage_pole).get_arc_direction(0)
# The angle that the segment deviates from the stage pole direction.
deviation_of_segment_direction_from_stage_pole = pygplates.Vector3D.angle_between(
segment_direction_at_midpoint, segment_to_stage_pole_direction)
# When comparing the deviation angle we need to consider the case where the two
# direction vectors are aligned but pointing in opposite directions.
if (deviation_of_segment_direction_from_stage_pole < isochron_segment_deviation_in_radians or
deviation_of_segment_direction_from_stage_pole > math.pi - isochron_segment_deviation_in_radians):
ridge_segments.append(segment)
else:
transform_segments.append(segment)
print 'Number ridge segments: %d' % len(ridge_segments)
print 'Number transform segments: %d' % len(transform_segments)
```

## Details

The rotations are loaded from a rotation file into a `pygplates.RotationModel`

.

```
rotation_model = pygplates.RotationModel('rotations.rot')
```

The isochron features are loaded into a `pygplates.FeatureCollection`

.

```
isochron_features = pygplates.FeatureCollection('isochrons.gpml')
```

The time period and conjugate plate IDs are obtained using `pygplates.Feature.get_valid_time()`

,
`pygplates.Feature.get_reconstruction_plate_id()`

and `pygplates.Feature.get_conjugate_plate_id()`

.

```
begin_time, end_time = isochron_feature.get_valid_time()
plate_id = isochron_feature.get_reconstruction_plate_id()
conjugate_plate_id = isochron_feature.get_conjugate_plate_id()
```

`begin_time`

of the isochron’s
plate `plate_id`

relative to its conjugate plate `conjugate_plate_id`

using
`pygplates.RotationModel.get_rotation()`

.`conjugate_plate_id`

. We could have
set it to zero and it shouldn’t change the result. We set it to the isochron’s conjugate plate just
in case there is no plate circuit path from plate zero to plate `conjugate_plate_id`

.```
stage_rotation = rotation_model.get_rotation(
begin_time + 1, plate_id, begin_time, conjugate_plate_id, conjugate_plate_id)
```

```
stage_pole, stage_angle_radians = stage_rotation.get_euler_pole_and_angle()
```

`pygplates.Feature`

the geometry will be in present day coordinates.…where \(P_{A}\) is the anchor plate.

Rearranging this gives us the rotation of moving plate \(P_{M}\) from present day to time \(t_{to}\):

Using the approach in Composing finite rotations we can write the transformation of a present day geometry on moving plate \(P_{M}\) to time \(t_{to}\) via the stage pole reference frame:

```
stage_pole_reference_frame = rotation_model.get_rotation(
begin_time, plate_id, 0, conjugate_plate_id, conjugate_plate_id)
stage_pole = stage_pole_reference_frame.get_inverse() * stage_pole
```

Next we iterate over the geometries of the isochron feature using `pygplates.Feature.get_geometries()`

.

Note

We are iterating over the un-rotated (or present day) geometries as noted above.

```
for isochron_geometry in isochron_feature.get_geometries():
```

We then iterate over the segments of the `polyline`

geometry
of the isochron using `pygplates.PolylineOnSphere.get_segments()`

.

```
for segment in isochron_geometry.get_segments():
```

`pygplates.PointOnSphere`

or a `pygplates.MultiPointOnSphere`

since those geometry types do not have segments.`pygplates.PolylineOnSphere(isochron_geometry).get_segments()`

instead of `isochron_geometry.get_segments()`

.A zero-length `segment`

has not direction so we ignore them.

```
if segment.is_zero_length():
continue
```

`pygplates.GreatCircleArc.get_arc_point()`

and
the segment direction (tangential to the globe) at the midpoint is found using
`pygplates.GreatCircleArc.get_arc_direction()`

```
segment_midpoint = segment.get_arc_point(0.5)
segment_direction_at_midpoint = segment.get_arc_direction(0.5)
```

Next we calculate a `3D vector`

from the segment mid-point towards
the stage pole by creating an `arc`

from the mid-point to the
stage pole and then getting the direction of the arc using `pygplates.GreatCircleArc.get_arc_direction()`

.

```
segment_to_stage_pole_direction = pygplates.GreatCircleArc(
segment_midpoint, stage_pole).get_arc_direction(0)
```

*from*the segment’s mid-point, but in different directions.

*radians*) between them is found using

`pygplates.Vector3D.angle_between()`

.```
deviation_of_segment_direction_from_stage_pole = pygplates.Vector3D.angle_between(
segment_direction_at_midpoint, segment_to_stage_pole_direction)
```

*isochron_segment_deviation_in_radians*parameter determines the maximum deviation angle at which at which the isochron segment switches from a ridge segment to a transform segment.

`math.pi - isochron_segment_deviation_in_radians`

is the threshold used when the isochron direction
is facing away from the stage pole (instead of towards it).```
if (deviation_of_segment_direction_from_stage_pole < isochron_segment_deviation_in_radians or
deviation_of_segment_direction_from_stage_pole > math.pi - isochron_segment_deviation_in_radians):
ridge_segments.append(segment)
else:
transform_segments.append(segment)
```