Refine plate polygons / lines below a predefined threshold for plotting

[brmather]

Some plotting artefacts occur when there are very few points that define a polygon or line. For example, mid-ocean ridge geometries in the plot below have very few points and do not plot correctly.

Maybe there is some way to do this in cartopy but I think we may need to densify these segments so that they follow (or closely approximate) the expected great circle geometry.

[tom-new]

The approach that was used to densify line segments/polygons (i.e. using shapely.segmentize will not correctly generate points which along great circles. shapely.segmentize operates on R², and so using it on geographic coordinates will create straight lines in lat–lon space, which in general will not create straight lines (great circles) on the surface of the sphere.

I needed to solve the same problem when I needed to render shapefiles from GPlates in ParaView, and at first I also tried to use shapely.segmentize. I outline the strategy I eventually used to correctly generate points along a great circle that connects an arbitrary pair of geographic coordinates below:

1. Convert the pair of geographic coordinates {a, b} to R³ such that they lie on the surface of the unit sphere.
2. Calculate the angle between them with the cross product θ = |a × b|.
3. Compare the angle to some threshold angle and calculate the number of intermediate vectors n that are required to sufficiently densify the line (might be zero, in which case we can skip all the next steps).
4. We now need basis vectors {x₀, x₁} in the plane of {a, b}. We can use a as x₀, and x₁ may be generated with the cross products
   x₂ = a × b / |a × b|, and
   x₁ = x₂ × x₀.
5. The n densifying vectors d_m may then be generated with the equation
   d_m = x₀cos(mθ / (n + 1)) + x₁sin(mθ / (n + 1)),
   where 1 ≤ m ≤ n.
6. Convert the densifying vectors back to geographic coordinates.

Below is the function I wrote in Python which performs steps 2 to 5 on a single pair of coordinates. It should be fairly easy to speed up by performing the calculations on shape (n,3) numpy.ndarray objects instead. The function also doesn't behave as wanted if the angle between a coordinate pair is reflex (i.e. < 180°). A priori this sounds like it could be annoying to solve, however it's also probably very rare.

import numpy as np

def densify_great_circle(array1: np.ndarray, array2: np.ndarray, max_angle: float=np.pi/180):
    '''
    Given two arrays of shape (3,) of points on the surface of a unit sphere, return a sequence of points along the great circle between them, with a separating angle of no more than `max_angle`.
    '''

    # get the angle between the two points
    angle = np.linalg.norm(np.cross(array1, array2))

    # exit if the angle between the points is already smaller than the maximum allowed angle
    if angle < max_angle:
        return

    # get the number of needed in-between vectors
    n = int(np.floor(angle / max_angle))

    # get the angle for the in-between vectors
    sub_angle = angle / (n + 1)

    # get basis vectors for plane on which the two vectors lie
    basis0 = array1 / np.linalg.norm(array1)
    basis2 = np.cross(array1, array2) / np.linalg.norm(np.cross(array1, array2))
    basis1 = np.cross(basis2, basis0)

    sub_vectors = []

    # create the in-between vectors
    for i in range(n):
        sub_vector = basis0 * np.cos((i + 1) * sub_angle) + basis1 * np.sin((i + 1) * sub_angle)
        sub_vector = sub_vector
        sub_vectors.append(sub_vector)

    return sub_vectors

[jcannon-gplates]

There's also pygplates.GreatCircleArc.to_tessellated() that I think should reproduce the above algorithm (ie, return uniformly spaced points along a great circle arc defined by two points).

And there's a version for polylines and polygons that repeats that process for each segment/arc of a polyline or polygon, but returns a tessellated polyline or polygon instead of a list of points (but that can be converted to a list of lat/lon points using polygeom.to_lat_lon_array())). Note that there can be non-uniform tessellation across great circle arcs - in the docs "The distance between adjacent points (in the tessellated polyline) will not be exactly uniform. This is because each segment in the original polyline is tessellated to the nearest integer number of points (that keeps that segment under the threshold) and hence each original segment will have a slightly different tessellation angle.".

By it may be easier to use your function. Especially if you want to change the details of the tessellation later - eg, uniform tessellation across an entire polyline or polygon - where you might specify an exact spacing angle (with no integer truncation) that leads to the points all being equidistant across the entire polyline/polygon (except for the last point). That's an option I will eventually add to pyGPlates.

[brmather]

The default tessellation is now set to 1 degree:

gplately/gplately/plot.py

Line 707 in bd3e4b8

tessellate_degrees = kwargs.pop("tessellate_degrees", 1)

These are not great circle distances, as far as I can see. The relevant code is in gplately.geometry.