GPlates / gplately

GPlately is a Python package to interrogate tectonic plate reconstructions.
https://gplates.github.io/gplately/
GNU General Public License v2.0
60 stars 13 forks source link

various polyline functions #251

Open michaelchin opened 2 months ago

michaelchin commented 2 months ago

These requests come from HS. DM wants these functions in gplately.

We need three functions:

1) Split Polylines into Segments Based on Azimuth Changes: This function breaks down polylines into segments based on changes in Azimuth. It converts lines into points, considering the distance between nodes (in degrees) and a threshold value (also in degrees) that determines the segmentation. If the azimuth change between a segment (which can consist of a line between two nodes or multiple connected nodes with minimal azimuth change) and the following subsegment (a line between two nodes) exceeds the threshold, the line is split.

2) Calculate Parameters for Line Segments: This function takes a set of points and polylines and computes various parameters for each point. These parameters include the distance from points to the nearest segment, as well as the azimuth, length, and velocity of the nearest segment.

3) Calculate Distance and Angle Between Segments: This function processes two sets of lines and a point dataset to calculate the distance between the closest segments to each point, as well as the angle between these segments. The distance between two lines can be measured in different ways, one of which is by averaging the distances between corresponding points on the two lines, including the starting point, midpoint, and endpoint.

michaelchin commented 2 months ago

some valuable comments from John https://github.com/GPlates/gplately/issues/255#issuecomment-2286795352

michaelchin commented 2 months ago

something we agreed https://github.com/GPlates/gplately/pull/258#issuecomment-2292840015

michaelchin commented 2 months ago

Hi @jcannon-gplates

I wonder if it is possible for you to help @Hojat-Shirmard with these polyline functionalities? The reason DM asked me to do this in the first place is because I had some old Python code which could do a small portion of this task. I have checked in my old code, and now I have to figure out how to do the rest. I reckon you are much more experienced with those geometry algorithms, such as calculating the angles, velocities, distance, etc. I guess you had done the similar things before with the subduction convergence stats and, hence, could do this task much faster and more correctly than I do.

Given enough time, I am happy to do these functionalities for @Hojat-Shirmard myself. But right now I have a few other things at hand to worry about, and, to make it worse, my teeth are acting up recently. So, I will very much appreciate it if you could help @Hojat-Shirmard and me out here. Thanks.

michaelchin commented 2 months ago

My old code to convert polyline into points has been checked in via https://github.com/GPlates/gplately/issues/255

jcannon-gplates commented 2 months ago

That's fine with me - yeah those functionalities are a little more up my alley. Let's check with Dietmar at the next meeting to be sure though (in light of his email just now). Hope your teeth get better!

michaelchin commented 2 months ago

The title of this issue is really bad. We should update the title once it is clear what to do. Maybe even split this issue into small ones.

jcannon-gplates commented 2 months ago

Hi @Hojat-Shirmard

I just want to clarify these functions some more:

Please note, I just realised that my definition of segment might be different than yours (in the following paragraphs). I'm defining it as a great circle arc between two consecutive points in a polyline. Whereas you might be defining it as more like a polyline that was split-off from a larger polyline? In any case, in the following I referred to a segment using my definition - just to be clear.

  1. Split Polylines into Segments Based on Azimuth Changes

Does this convert polylines into smaller polylines ? Or does it convert into a list of point runs (ie, a list of a list of points) where each point run is a continuous part of the original polyline that has been split off?

And it also involves tessellation, right? As in, a new point is created every threshold distance along the polyline, for each point run (starting at the first point). And what about the last point in a single point run - ie, the last point on the original polyline segment before a split point (which is also the first point on the next point run) - should it be included in the point run?

And a new point run is started whenever a segment (defined as a single great circle arc between two original points in the original polyline) changes direction enough relative to the previous segment (in original polyline)? I assume by azimuth you're just referring to the direction (which is not necessarily relative to North, ie, is just relative to the previous segment).

And when you say If the azimuth change between a segment (which can consist of a line between two nodes or multiple connected nodes with minimal azimuth change - why is the part in bold needed? Can it just be determined by the two original points of the segment in question (rather than also looking at previous segments)? As in, you just look at the direction of the following segment relative to the current segment (ignoring previous segments) to determine when to split. By the way, by nodes I assume you're referring to the original points in the original polyline, is that right? As in, a segment in the original polyline is between two nodes.

  1. Calculate Parameters for Line Segments

So the distance from points to the nearest segment is a distance, for each point, to the nearest polyline? By segment it sounds like you're referring to one of the input polylines, but I assume you just mean a great circle arc of a polyline because you also mention azimuth and length of nearest segment.

The velocity is a little trickier because this is the only request that does not purely involve geometry. Ie, it depends on where the geometry came from (eg, a line reconstructed by plate ID, or by half-stage rotations, or maybe it's a resolved topology). So I think we should not include that here. There are plans to define velocities along various reconstructed/resolved entities in pyGPlates - so we should probably wait for that and then build functions in GPlately from that (at a higher code level than pure geometry, where there's some context about reconstruction). Or the user can calculate velocities manually in the meantime.

  1. Calculate Distance and Angle Between Segments

I assume this takes each point in a point list and finds the nearest line in line set 1 (say at point A on that nearest line) and the nearest line in line set 2 (say at point B on that nearest line) and then calculates the distance between points A and B? You mention different ways to calculate the distance, but I'd imagine just calculating between A and B is sufficient - pyGPlates can calculate the closest distance from a point to a line (where point A on the line might be at a vertex of the line or even between vertices) - see pygplates.GeometryOnSphere.distance(). You also mention the angle between these segments - but I'm not sure what angle that is?

Hojat-Shirmard commented 2 months ago

Hi @jcannon-gplates,

We need to raise these questions with Ehsan (e-farahbakhsh). Can you please add him here?

jcannon-gplates commented 2 months ago

Hi @jcannon-gplates,

We need to raise these questions with Ehsan (e-farahbakhsh). Can you please add him here?

No problem, I've sent out an invite.

jcannon-gplates commented 2 months ago

Hi @e-farahbakhsh, Hojat requested your input on my questions above. No hurry, just letting you know.

Hojat-Shirmard commented 2 months ago

Hi @jcannon-gplates, I successfully wrote the first function (1-Split Polylines into Segments Based on Azimuth Changes)!

jcannon-gplates commented 1 month ago

I successfully wrote the first function (1-Split Polylines into Segments Based on Azimuth Changes)!

Great! Could you perhaps paste it in a comment here? Only for reference, it might help answer some of the questions above. And ultimately we'd implement a version in gplately.

Hojat-Shirmard commented 1 month ago

Yes, absolutely. I'm currently working on the other functions, and once I've completed them all, I'll share everything with you in a single file for implementation in gplately.

Hojat-Shirmard commented 1 month ago

Here is the summary of all three functions that I prepared in a single notebook.

Calculate Spatial Relationships

Split polylines into segments based on azimuth changes. Calculate parameters for line segments and assign values to nearest points, including azimuth, length, and the nearest distance to points (from start to endpoints of segments). Calculate the distance and angle between two segments of lines and assign values to the nearest points.
Inputs are two lines (here are the craton and rift lines) and one point (here is the deposit location) in ".shp" format.
At the end, we will extract 8 key parameters:

1- craton_dis: The minimum distance from craton segments (line 1) to the nearest deposit (point).
2- craton_az: The azimuth of craton segments (line 1), assigned to the nearest deposit (point).
3- craton_len: The length of craton segments (line 1), assigned to the nearest deposit (point).
4- rift_dis: The minimum distance from rift segments (line 2) to the nearest deposit (point).
5- rift_az: The azimuth of rift segments (line 2), assigned to the nearest deposit (point).
6- rift_len: The length of rift segments (line 2), assigned to the nearest deposit (point).
7- seg_dis: The distance between craton and rift segments (line 1 and line 2), assigned to the nearest deposit (point).
8- seg_angle: The smallest angle between craton and rift segments (line 1 and line 2), assigned to the nearest deposit (point).

These parameters provide valuable spatial relationships between geological features (lines) and deposits (points).

The final output will be: "merged_segmented_rift" (line 1), "merged_segmented_craton" (line 2), and "Deposits_SpatialRelationships" (point).
"Deposits_SpatialRelationships" (point) will carry all 8 extracted parameters