How to create a small circle features in pygplates

Q&A
1

Hello everyone: may I ask whether it is possible to create small circle features in pygplates similarly like how one can create small circle features in GPlates?

Thank you for your time and your help.

2

Hi Lavie,

You can create a small circle feature in pyGPlates using the following code…

import pygplates

# Small circle parameters.
small_circle_latitude = 0.0
small_circle_longitude = 0.0
small_circle_angular_radius_degrees = 15.0
small_circle_reconstruction_plate_id = 701
small_circle_name = 'Small Circle'
small_circle_begin_time = 100.0
small_circle_end_time = 0.0

# Create empty small circle feature.
small_circle_feature_type = pygplates.FeatureType.create_gpml('SmallCircle')
small_circle_feature = pygplates.Feature(small_circle_feature_type)

# Small circle centre point.
small_circle_centre = pygplates.PointOnSphere(small_circle_latitude, small_circle_longitude)

# Add parameters to small circle feature.
small_circle_feature.set_geometry(small_circle_centre)
small_circle_feature.set_reconstruction_plate_id(small_circle_reconstruction_plate_id)
small_circle_feature.set_name(small_circle_name)
small_circle_feature.set_valid_time(small_circle_begin_time, small_circle_end_time)
# Setting the angular radius parameter is more complicated:
# 1) You need to explicitly specify the property name 'gpml:angularRadius'.
#    This is because there's no common function like 'small_circle_feature.set_angular_radius()'.
# 2) You need to turn off validation with the GPlates Geological Information Model (GPGIM).
#    This is because the GPGIM expects the property to be a GpmlMeasure
#    (which we don't have in pyGPlates). Instead we create essentially the same thing
#    (an XsDouble property) using 'small_circle_feature.set_double()'.
#    But we have to tell the GPGIM to disregard this.
small_circle_feature.set_double(
        pygplates.PropertyName.create_gpml('angularRadius'),
        small_circle_angular_radius_degrees,
        pygplates.VerifyInformationModel.no)

# Save the small circle feature to a GPML file.
pygplates.FeatureCollection(small_circle_feature).write('small_circle_pygplates.gpml')

…note that setting the angular radius is more complicated than normal (see the comment in the code) but the resulting GPML file still loads into GPlates and displays the small circle correctly.

The code essentially involves creating an empty small circle feature and then populating it with the small circle parameters like the centre point and angular radius, as well as the more common properties like reconstruction plate ID, valid time period and name.

A slight alternative is to use pygplates.create_reconstructable_feature…

import pygplates

# Small circle parameters.
small_circle_latitude = 0.0
small_circle_longitude = 0.0
small_circle_angular_radius_degrees = 15.0
small_circle_reconstruction_plate_id = 701
small_circle_name = 'Small Circle'
small_circle_begin_time = 100.0
small_circle_end_time = 0.0

# Small circle feature type.
small_circle_feature_type = pygplates.FeatureType.create_gpml('SmallCircle')

# Small circle centre point.
small_circle_centre = pygplates.PointOnSphere(small_circle_latitude, small_circle_longitude)

# Create a feature with the common parameters in one call.
small_circle_feature = pygplates.Feature.create_reconstructable_feature(
        small_circle_feature_type,
        small_circle_centre,
        name=small_circle_name,
        valid_time=(small_circle_begin_time, small_circle_end_time),
        reconstruction_plate_id=small_circle_reconstruction_plate_id)
# Add the remaining (more complicated) angular radius parameter.
#
# Setting the angular radius parameter is more complicated:
# 1) You need to explicitly specify the property name 'gpml:angularRadius'.
#    This is because there's no common function like 'small_circle_feature.set_angular_radius()'.
# 2) You need to turn off validation with the GPlates Geological Information Model (GPGIM).
#    This is because the GPGIM expects the property to be a GpmlMeasure
#    (which we don't have in pyGPlates). Instead we create essentially the same thing
#    (an XsDouble property) using 'small_circle_feature.set_double()'.
#    But we have to tell the GPGIM to disregard this.
small_circle_feature.set_double(
        pygplates.PropertyName.create_gpml('angularRadius'),
        small_circle_angular_radius_degrees,
        pygplates.VerifyInformationModel.no)

# Save the small circle feature to a GPML file.
pygplates.FeatureCollection(small_circle_feature).write('small_circle_pygplates.gpml')

Regards,
John

3

Dear John:

Whoa, that is amazing :slight_smile:

I very much appreciate for your time and your help. It is very helpful. Thank you so much.

4

Dear John:

Is it possible to get the coordinate of each point forming the small circle feature (like a PolylineOnSphere.to_lat_lon_list())? When I tried ‘get_geometry()’ for the small circle feature, I got the center of the circle. When I tried ‘get_property()’ for the small circle feature, I got the angular radius.

Many thanks.

5

Hi Lavie,

Not directly. As you noticed, small circle features store a point and a radius from which GPlates internally generates a tessellated small circle geometry at render time (in the globe/map views). So the actual small circle geometry is not directly available in the small circle feature (eg, in the GPML you only see the centre point and angular radius). This is also the case if you subsequently export it in GPlates.

However this pyGPlates code does essentially what GPlates does internally. So you can use this to generate a small circle geometry (as a polygon) with however many points you want (see the num_points_in_small_circle parameter). You can just add this code to the code I posted above…


import math

def get_any_perpendicular(point_on_sphere):
    """
    Find any point perpendicular to 'point_on_sphere'.
    In other words any point on the great circle that is 90 degrees away from 'point_on_sphere'.
    """
    # The vector cross product of 'point_on_sphere' with any point that is not coincident with it will point
    # in a direction that is perpendicular to 'point_on_sphere'.
    unnormalised_perp = pygplates.Vector3D.cross(pygplates.Vector3D.x_axis, point_on_sphere.to_xyz())
    if unnormalised_perp.is_zero_magnitude():
        # Oops, seems the 'point_on_sphere' is along the x-axis, so pick the y-axis instead.
        unnormalised_perp = pygplates.Vector3D.cross(pygplates.Vector3D.y_axis, point_on_sphere.to_xyz())
    
    # Points on the sphere are unit length (normalisation will make it unit length).
    return pygplates.PointOnSphere(unnormalised_perp.to_normalised().to_xyz())


def tessellate_small_circle(
        centre_point_on_sphere,
        angular_radius_degrees,
        num_points):
    """
    Generate 'num_points' points (pygplates.PointOnSphere) for the small circle with specified centre and angular radius.
    """
    # Get any point on the small circle.
    # This will be our starting point for generating the remaining small circle points.
    #
    # The start point is obtained by taking any rotation axis perpendicular to the small circle centre and
    # rotating the small circle centre, using that rotation axis, by the small circle angular radius.
    start_point_rotation_axis = get_any_perpendicular(centre_point_on_sphere)
    start_point_rotation = pygplates.FiniteRotation(
        start_point_rotation_axis,
        math.radians(angular_radius_degrees))
    start_point = start_point_rotation * centre_point_on_sphere
    
    # Get the rotation that generates the next tessellated point in the small circle (from the previous point).
    tessellate_angle_radians = 2 * math.pi / num_points
    tessellate_rotation = pygplates.FiniteRotation(
        centre_point_on_sphere,
        tessellate_angle_radians)
    
    small_circle_points = []
    
    # Starting  with the start point (on the small circle) generate all 'num_points' points on the small circle.
    small_circle_point = start_point
    for n in range(num_points):
        small_circle_point = tessellate_rotation * small_circle_point
        small_circle_points.append(small_circle_point)
    
    return small_circle_points


num_points_in_small_circle = 100

#print(get_any_perpendicular(small_circle_centre).to_lat_lon())
small_circle_points = tessellate_small_circle(
        small_circle_centre,
        small_circle_angular_radius_degrees,
        num_points_in_small_circle)

# Create a small circle feature, but of type 'UnclassifiedFeature' since type 'SmallCircle'
# is expecting a centre *point* and angular radius. But we're giving it a polygon
# (tessellated small circle). Note that the polygon is just to join the dots of our
# tessellated small circle (each segment/arc of the polygon is actually a *great* circle arc,
# but when you have enough points you don't notice this).
small_circle_polygon = pygplates.PolygonOnSphere(small_circle_points)
small_circle_tessellated_feature = pygplates.Feature.create_reconstructable_feature(
        pygplates.FeatureType.gpml_unclassified_feature,
        small_circle_polygon,
        name=small_circle_name,
        valid_time=(small_circle_begin_time, small_circle_end_time),
        reconstruction_plate_id=small_circle_reconstruction_plate_id)

# Save the tessellated small circle feature to a GPML file.
pygplates.FeatureCollection(small_circle_tessellated_feature).write('small_circle_tessellated_pygplates.gpml')

If you load the generated GPML file into GPlates you should see something like this (where I use the Move Vertex tool to highlight the points)…

image
image654×604 8.99 KB

Regards,
John

6

Dear John:

Your answer is packed with details and very helpful; I really do appreciate for your help.

Best,