Module gplately.ptt.utils.points_in_polygons

Efficient point-in-polygon testing when there are many relatively uniformly spaced points to be tested against polygons.

For example, to find the plate ID of the polygon containing each point in a sequence of points:

import points_in_polygons

# A list of 'pygplates.PointOnSphere' points.
points = [...]

# Some polygon features (eg, coastlines).
polygon_feature_collection = pygplates.FeatureCollection('polygons.gpml')

# Extract the polygons from the features.
polygons = []
polygon_features = []
for polygon_feature in polygon_feature_collection:
    polygons.append(polygon_feature.get_geometry())
    polygon_features.append(polygon_feature)

#
# Assuming the polygons are *non-overlapping*, find the single polygon (feature) containing each point.
#
polygon_features_containing_points = points_in_polygons.find_polygons(points, polygons, polygon_features)

# Assign a plate ID to each point (or 0 if point outside all polygons).
plate_ids = []
for polygon_feature in polygon_features_containing_points:
    if polygon_feature is not None:
        plate_id = polygon_feature.get_reconstruction_plate_id()
    else:
        plate_id = 0

    plate_ids.append(plate_id)

#
# Assuming the polygons are *overlapping*, find all polygons (features) containing each point.
#
polygon_features_containing_points = points_in_polygons.find_polygons(points, polygons, polygon_features, all_polygons=True)

# Assign multiple points to each plate ID.
plate_id_to_points_mapping = {} # Each plate ID has a list of points.
for point_index, polygon_feature_list in enumerate(polygon_features_containing_points):
    if polygon_feature_list:
        for polygon_feature in polygon_feature_list
            plate_id = polygon_feature.get_reconstruction_plate_id()
            points_with_plate_id = plate_id_to_points_mapping.setdefault(plate_id, [])
            points_with_plate_id.append(points[point_index])

Functions

def find_polygons(points, polygons, polygon_proxies=None, all_polygons=False, subdivision_depth=4)

Efficient point-in-polygon testing when there are many relatively uniformly spaced points to be tested against polygons.

Parameters

points : a sequence of 'pygplates.PointOnSphere'

a sequence of points

polygons : a sequence of 'pygplates.PolygonOnSphere'

a sequence of polygons

polygon_proxies : Optional sequence of objects associated with 'polygons'

If not specified then the proxies default to the polygons themselves. These can be any object (such as the 'pygplates.Feature' that the polygon came from).

all_polygons : bool

Whether to find all polygons containing each point or just the first one encountered. Set to True if polygons overlap each other, otherwise set to False (for non-overlapping polygons). Defaults to False (non-overlapping polygons).

subdivision_depth : number

The depth of the lat/lon quad tree used to speed up point-in-polygon queries. The lat/lon width of a leaf quad tree node is (90 / (2^subdivision_depth)) degrees. Generally the denser the 'points' the larger the depth should be. Setting this value too high causes unnecessary time to be spent generating a deep quad tree. Setting this value too low reduces the culling efficiency of the quad tree. However a value of 4 seems to work quite well for a uniform lat/lon spacing of 'points' of 1 degree and below without the cost of generating a deep quad tree. So most of the time the subdivision depth can be left at its default value.

Returns

A list of polygon proxies associated with 'points'

The length of the returned list matches the length of 'points'. For each point in 'points', if the point is contained by a polygon then that polygon's proxy is stored (otherwise None is stored) at the same index (as the point) in the returned list. If 'all_polygons' is False then each item in returned list is a single polygon proxy (or a single None). If 'all_polygons' is True then each item in returned list is a list of polygon proxies (or a single None).

Raises ValueError if the lengths of 'polygons' and 'polygon_proxies' (if specified) do not match.

def find_polygons_using_points_spatial_tree(points, spatial_tree_of_points, polygons, polygon_proxies=None, all_polygons=False)

Same as 'find_polygons()' except 'spatial_tree_of_points' is a 'points_spatial_tree.PointsSpatialTree' of 'points'.

This is useful when re-using a single 'points_spatial_tree.PointsSpatialTree'. For example, when using it both for point-in-polygon queries and minimum distance queries.

Note that 'spatial_tree_of_points' should have been built from 'points' since it contains indices into the 'points' sequence.