Find nearest reconstructed feature to a point

This example iterates over a collection of reconstructed features and finds the feature that is nearest to a point.

Sample code

import pygplates


# Load one or more rotation files into a rotation model.
rotation_model = pygplates.RotationModel('rotations.rot')

# Load some features.
features = pygplates.FeatureCollection('features.gpml')

# Reconstruct features to 10Ma.
reconstruction_time = 10

# All features have their distance calculated relative to this point.
point_latitude = 0
point_longitude = 0
point = pygplates.PointOnSphere(point_latitude, point_longitude)

# Reconstruct the features to 10Ma.
reconstructed_features = []
pygplates.reconstruct(features, rotation_model, reconstructed_features, reconstruction_time, group_with_feature=True)

# The minimum distance to all features and the nearest feature.
min_distance_to_all_features = None
nearest_feature = None

# Iterate over all reconstructed features.
for feature, feature_reconstructed_geometries in reconstructed_features:

    # Iterate over all reconstructed geometries of the current feature.
    for feature_reconstructed_geometry in feature_reconstructed_geometries:

        # Get the minimum distance from point to the current reconstructed geometry.
        min_distance_to_feature = pygplates.GeometryOnSphere.distance(
                point,
                feature_reconstructed_geometry.get_reconstructed_geometry(),
                min_distance_to_all_features)

        # If the current geometry is nearer than all previous geometries then
        # its associated feature is the nearest feature so far.
        if min_distance_to_feature is not None:
            min_distance_to_all_features = min_distance_to_feature
            nearest_feature = feature

if nearest_feature is not None:
    print "The nearest feature, to point %s, has feature ID %s and a minimum distance of %fkms" % (
            point.to_lat_lon(),
            nearest_feature.get_feature_id(),
            min_distance_to_all_features * pygplates.Earth.mean_radius_in_kms)

Details

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

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

The reconstructable features are loaded into a pygplates.FeatureCollection.

features = pygplates.FeatureCollection('features.gpml')

The features are reconstructed to their 10Ma positions.

reconstruction_time = 10
The test point has zero latitude and zero longitude.
All features are tested to see which one is closest to this point.
point_latitude = 0
point_longitude = 0
point = pygplates.PointOnSphere(point_latitude, point_longitude)
All features are reconstructed to 10Ma using pygplates.reconstruct().
We specify a list for reconstructed_features instead of a filename.
We also set the output parameter group_with_feature to True (it defaults to False) so that our reconstructed feature geometries are grouped with their feature. This isn’t strictly necessary in this particular example though.
reconstructed_features = []
pygplates.reconstruct(features, rotation_model, reconstructed_features, reconstruction_time, group_with_feature=True)
Initially we don’t have a minimum distance of the point to all features.
This value is also used as the threshold to the distance function and initially this will be None which means no threshold.
min_distance_to_all_features = None
nearest_feature = None

Each item in the reconstructed_features list is a tuple containing a feature and its associated reconstructed geometries.

for feature, feature_reconstructed_geometries in reconstructed_features:

A feature can have more than one geometry and hence will have more than one reconstructed geometry.

for feature_reconstructed_geometry in feature_reconstructed_geometries:
Calculate the minimum distance from the point to a reconstructed feature geometry using pygplates.GeometryOnSphere.distance().
min_distance_to_all_features is specified as the distance threshold since we’re only interested in geometries that are nearer than the closest geometry encountered so far.
min_distance_to_feature = pygplates.GeometryOnSphere.distance(
        point,
        feature_reconstructed_geometry.get_reconstructed_geometry(),
        min_distance_to_all_features)
If None was returned then the distance was greater than min_distance_to_all_features.
So a valid returned value means the current geometry is the nearest geometry encountered so far.
In this case we record the nearest feature and the new minimum distance.
if min_distance_to_feature is not None:
    min_distance_to_all_features = min_distance_to_feature
    nearest_feature = feature

Once we’ve tested all features (if any were in the file) we print out the nearest feature and its (minimum) distance to the point.

if nearest_feature is not None:
    print "The nearest feature, to point %s, has feature ID %s and a minimum distance of %fkms" % (
            point.to_lat_lon(),
            nearest_feature.get_feature_id(),
            min_distance_to_all_features * pygplates.Earth.mean_radius_in_kms)

Advanced

If we also want to find the nearest position on the nearest feature then we can specify return_closest_positions as True when calculating the geometry distance:

nearest_point_on_nearest_feature = None

...

min_distance_to_feature_and_closest_positions = pygplates.GeometryOnSphere.distance(
        point,
        feature_reconstructed_geometry.get_reconstructed_geometry(),
        min_distance_to_all_features,
        return_closest_positions=True)
When return_closest_positions is True the distance function returns a 3-tuple of (distance, closest point on geometry1, closest point on geometry2).
If the geometry distance is below the threshold then we can extract that tuple into its individual elements.
if min_distance_to_feature_and_closest_positions is not None:
    min_distance_to_all_features, _, nearest_point_on_nearest_feature = min_distance_to_feature_and_closest_positions
    nearest_feature = feature
…where the _ is just a placeholder to ignore the second element of the tuple returned by the distance function.
We ignore it because it’s the closest distance to the test point which is just the test point itself. And we are only interested in the closest point on the feature’s geometry.