unify the geometries of COBs in Muller2016

[michaelchin]

Currently, the geometries of COBs in Muller2016 are a mixture of polygons and polylines. The mixture is causing problem when being plotted in GPlately. We need to modify the COBs and make the geometries all polygons(or polylines?).

Thanks to @nickywright for providing us with a code snippet to convert polylines to polygons.

import os.path
import pygplates


input_filenames = [
    'Global_EarthByte_GPlates_PresentDay_StaticPlatePolygons_2019_v1.shp'
]

output_filename_append = '_polygon'

for input_filename in input_filenames:
    # Convert polylines to polygons in the features.
    feature_collection = pygplates.FeatureCollection(input_filename)
    for feature in feature_collection:
        # polygons = [pygplates.PolygonOnSphere(geometry) for geometry in feature.get_geometries()]
        # changed feature get statement to 'all_geometries' - since feature has issues getting different
        # types of geometries
        for geometry in feature.get_all_geometries():
            polygons = pygplates.PolygonOnSphere(geometry)
            feature.set_geometry(polygons)
    
    # Get the output (polygon) filename from the input (polyline) filename.
    input_filename_root, input_filename_ext = os.path.splitext(input_filename)
    output_filename = ''.join((input_filename_root, output_filename_append, input_filename_ext))
    
    # Write the modified feature collection to output file.
    pygplates.FeatureCollection(feature_collection).write(output_filename)

Thanks to @jcannon-gplates for providing valuable suggestion.

basically a pygplates.PolygonOnSphere can take a polyline as a single argument and it’ll convert it (see the second init() listed here). And similarly for converting a polygon (exterior) to a polyline. So if you can get the polyline/polygon into pygplates format then converting is straightforward.

[jcannon-gplates]

I see that it’s been modified to use get_all_geometries () instead of get_geometries() – the only issue might be if have more than one geometry in a single feature - it'll set the last one only. But that should be fine since usually only one geometry per feature anyway.

If it became an issue you could instead have something like:

polygons = [pygplates.PolygonOnSphere(geometry) for geometry in feature.get_all_geometries()]
feature.set_geometry(polygons)

[brmather]

Good pick up, @jcannon-gplates!

[brmather]

Has there been any progress made on this?

[michaelchin]

Has there been any progress made on this?

Can we find a RA to do this?

[michaelchin]

The COBs in GeoData 2.5 are polylines. Shall we follow it as example and convert all COB polygons in Muller2016 to polylines?

Alternatively, we can create two COB files for Muller2016, one with polygons and the other with polylines.

Hi @bsim8882, Just checking if you can see the email notification from github.com. If you see this, please let me know. We will use github.com extensively to discuss issues. It is important that you can see the notifications.

[michaelchin]

Let's follow John's suggestion https://github.com/GPlates/gplately/pull/302#issuecomment-2499613060

I think we should allow users to do both way, polyline->ploygon and polygon->polyline.

[nickywright]

@michaelchin, Just my 2c: need to check what the COBs actually look like, because maybe the COBs from the GeoData should stay as polylines. Converting these to polygons is going to create lots of issues because these aren't in anyway closed like this (the file from the COB folder in 2.5):

I use the script above to convert the agegridding-masking COBs to polygons for plotting/analysis. I don't know if GPlately has these at all, but for example, for Müller et al. 2019 it is the file on the USyd-github in EarthBytePlateMotionModel-ARCHIVE/Muller_++_2019_PlateMotionModel_Tectonics/StaticGeometries/AgeGridInput/Global_EarthByte_GeeK07_COB_Terranes_2019_v2.gpml and looks like this (and I think there is a couple of pesky polylines in it amongst all the polygons):

[michaelchin]

Thanks @nickywright

Yes, you are right. We should not change the COBs in GeoData. We will mainly focus on fixing the COBs in Muller2016 in this Issue.

For model Muller2019, we are using a file named Global_EarthByte_GeeK07_COBLineSegments_2019_v1.gpmlz and it looks similar to the one in GeoData. I wonder if there is such a file in model Muller2016 as well.

Anyway, it is no harm to add new functions to convert between polylines and polygons in GPlately. And Bianca can use this Issue to get started with GPlates development work.

[nickywright]

I agree that adding the polylines-->polygons function would be super handy!

[jcannon-gplates]

Yes I think converting polygons to polylines is fine. But converting polylines to polygons has the problem Nicky mentioned (ie, polylines aren't closed).

One idea for conversion script is to detect if the polylines are closed and only convert if so. I vaguely recall some old datasets that ended up being closed polylines instead of polygons. These would be good for converting to polygons.

Some random issues that may arise when converting an entire dataset to a single geometry type as required for #214:

• maybe conversion has to be all or nothing (since don't want to end up with a mix of geometry types)
• what if points/multi-points are also in the dataset (they probably shouldn't be converted to polylines/polygons)

[michaelchin]

Yes, I agree that we need a comprehensive plan down the road. For this Issue, I would like to focus on fixing the COBs for Muller2016 model only. We can create new github Issues to track the various ideas. But I would like to make this one simple and clear. I don't want Bianca feel overwhelmed at this stage.

[jcannon-gplates]

No problem 👍. Yes, let's wait and see how this pans out. And then later we can create a new issue if something more comprehensive is needed.

[michaelchin]

@jcannon-gplates I might have done something wrong or misunderstood the feature.set_geometry(polygons) function. See the debug code below. It seemed the set_geometry() just added the converted geometry and did not remove the old geometries. Did I do something wrong?

https://github.com/GPlates/gplately/blob/86559303f9bae60dda9c0458a879c53b723cbe3d/gplately/utils/convert_geometries.py#L36-L41

[michaelchin]

@bsim8882 Could you please try to follow the instructions in the source code? Thanks.

https://github.com/GPlates/gplately/blob/86559303f9bae60dda9c0458a879c53b723cbe3d/tests-dir/unittest/test_convert_geometries.py#L27-L44

[jcannon-gplates]

@jcannon-gplates I might have done something wrong or misunderstood the feature.set_geometry(polygons) function. See the debug code below. It seemed the set_geometry() just added the converted geometry and did not remove the old geometries. Did I do something wrong?

I think your understanding is correct.

I suspect what is happening is the property name of the original polyline is not the default geometry property name.

And so set_geometry() is setting the converted polygon into a new property (leaving the original polyline property untouched). If they had the same property name then the original polyline property would have been replaced. So I think this issue is caused by the original polyline being in an unexpected property name (according to the GPGIM). Note that set_geometry() has a second argument that is the property name (if not specified then assumes the default geometry property name).

This is why I noted to myself in the PR here that...

Also should retain property name of each geometry property. For inspiration, reference this.

...because I think using get_all_geometries() in this manner is quite risky. It just scoops up all the geometries regardless of their name or type. For example, it could grab some point geometry that represents a centroid of the feature (in some weird feature type), but that feature could also have a regular polyline property under a different property name. And then we'd be setting the point and polyline under the same property name (essentially corrupting that feature in terms of the GPGIM).

It's actually a little bit tricky to do this properly. Perhaps something like this, which converts any polylines to polygons (in a single feature), retaining their property names and leaving non-polylines unchanged:

# First step: gather all geometry property names (ignore other properties).
geometry_property_names = []
for property in feature:
    # If it's a geometry property then add the property name to the list.
    if property.get_value().get_geometry():
        geometry_property_names.append(property.get_name())

# Second step: convert polylines to polygons (don't convert other geometry types).
for geometry_property_name in geometry_property_names:
    # Get all geometries with the current property name.
    #
    # NOTE: I'm specifying the property name rather than relying on the default.
    geometries_with_property_name = feature.get_geometries(geometry_property_name)

    # Among the current geometries, convert any polylines to polygons (leave the other types alone).
    for geometry_index, geometry in enumerate(geometries_with_property_name):
        if isinstance(geometry, pygplates.PolylineOnSphere):
            # As a bonus, only convert polyline if it is closed (ie, its first and last vertices are equal).
            # Otherwise converting to a polygon will be dodgy (as Nicky pointed out).
            if geometry[0] == geometry[-1]:
                polygon = pygplates.PolygonOnSphere(geometry)
                geometries_with_property_name[geometry_index] = polygon

    # Set all geometries with the current property name.
    #
    # NOTE: I'm specifying the property name rather than relying on the default.
    feature.set_geometry(geometries_with_property_name, geometry_property_name)

...where the first step gathers all the property names associated with geometry properties (ie, ignores plate ID properties, etc). And the second step only converts polylines to polygons (and only if the polylines are closed). Also, it doesn't change the property names inadvertently. And it retains any original points, multi-points and polygons.

A similar thing could be done for converting the other way - from polygons to polylines (except without checking that polygons are closed).

A final note, this should pretty much avoid any GPGIM issues, but there's still a small chance it could happen if a particular geometry property name only supports a subset of geometry types (eg, might only support polygons, not polylines, and so can't convert polygon to polyline). But I don't think any of these properties actually exist in the GPGIM currently (for example, the gpml:centerLineOf property supports all geometry types). In any case, if there was a problem then the verify_information_model argument of set_geometry() could be turned off (it's on by default).

[michaelchin]

hi @bsim8882

commit 926b315c4e111cd4f4d207fd2b3c27c6b7cfa58f looks great. Just one tiny thing, normally we don't check in data/log/generated files. Ideally, we only check in source code/configuration files(basically files with a good reason to track).

And do you want to try do this function?

https://github.com/GPlates/gplately/blob/0a72e3f1c49fcecb8467805ad76084a0f3da2a04/gplately/utils/convert_geometries.py#L83-L103

@jcannon-gplates Do you have better ideas of converting polygons to polylines?

[jcannon-gplates]

# step 2: get the exterior ring and interior rings from the polygon and convert them to PolylineOnSphere
#         John may have better ideas on this?
#         see the links below for relavant pygplates doc

There's a tricky gotcha here.

If you convert using:

polyline = pygplates.PolylineOnSphere(polygon)

...then it'll convert the exterior ring only (ignoring any interior rings).

But it'll also make sure the converted polyline is closed (because it knows it's a polygon).

Alternatively, if you do this:

polyline = pygplates.PolylineOnSphere(polygon.get_exterior_ring_points())

...then it won't ensure the converted polyline is closed (because it doesn't know the points are from a polygon).

For example, if you convert a triangle polygon (with vertices A, B, C) to a polyline then:

• in the first case you'll get 4 vertices A->B->C->A (ie, the polyline with have the 3 segments of the triangle).
• in the second case you'll get 3 vertices A->B->C only (ie, only two segments, you'll be missing the segment C->A).

So I'd recommend the first case if we're only converting the exterior ring of a polygon to a polyline.

[jcannon-gplates]

Then there's the question of whether to also make polylines from interior rings?

To be honest, I'm not sure. Maybe @nickywright has an idea on whether there are datasets where, say, a polygon with 2 interior rings should be converted to 3 separate polylines (one for the exterior ring and two for the interior rings)?

If so, then I can add to pyGPlates the ability to do:

polylines = pygplates.PolylineOnSphere(polygon, convert_interior_rings=True)

...where specifying convert_interior_rings=True will return a list of polylines (versus the default convert_interior_rings=False remaining unchanged; ie, returning a single polyline, not a list).

Update: Actually it'd probably be more like:

polylines = pygplates.PolylineOnSphere.convert_from(polygon, convert_interior_rings=True)

...since pygplates.PolylineOnSphere(...) should always return a single polylline.

[michaelchin]

Hi @bsim8882

let's do what John said above. It is simple and clean.

polyline = pygplates.PolylineOnSphere(polygon)

[brmather]

Does converting polygons to polylines impact on the masking of age grids?

[michaelchin]

Does converting polygons to polylines impact on the masking of age grids?

I am unsure about this. @jcannon-gplates do you know anything about this?

[jcannon-gplates]

For regular masking (ie, not continent contouring)...

It looks like polylines are converted to polygons in ensure_polygon_geometry() here: https://github.com/GPlates/gplately/blob/ece4c2d059266f9475daa658b600a5733dbc7f21/gplately/utils/seafloor_grid_utils.py#L117

...so that's good. However, ensure_polygon_geometry() is not always used for masking, as far as I can tell. For example...

I'm not sure what rasterio.features.rasterize will do with a polyline? It probably won't fill it. https://github.com/GPlates/gplately/blob/ece4c2d059266f9475daa658b600a5733dbc7f21/gplately/oceans.py#L802-L803

...it may be that it should instead be called with self.continental_polygons (rather than self._PlotTopologies_object.continents)? ...which is generated by the aforementioned ensure_polygon_geometry() here: https://github.com/GPlates/gplately/blob/ece4c2d059266f9475daa658b600a5733dbc7f21/gplately/oceans.py#L296-L298

[jcannon-gplates]

For continent contouring...

Luckily it looks like polylines are converted to polygons before the contouring begins: https://github.com/GPlates/gplately/blob/ece4c2d059266f9475daa658b600a5733dbc7f21/gplately/ptt/continent_contours.py#L641-L653

As noted, it only works if the polylines are closed. But if we've converted polygons to polylines beforehand then they will be closed polylines.

[jcannon-gplates]

Does converting polygons to polylines impact on the masking of age grids?

So overall it might be ok with the above-mentioned fixup for regular (non-contoured) masking.

The only other issue is polygons with interior holes would not have the holes contribute to the mask.