Let gmtspatial -W extend segments from one of both ends #6255

PaulWessel · 2022-01-23T03:15:11Z

This PR follows the comments on the forum post and adds a new -W option that is used to extend existing lines in one or both directions by a given amount(s). A test has been added to show how it works for geographic and Cartesian settings, and the documentation has been updated. Here is the test output:

This PR follows the comments on the forum post and adds a new -W option that is used to extend existing lines in one or both directions by a given amount. A test has been added and the documentation has been updated.

Esteban82

Many thanks Paul!!

Esteban82 · 2022-01-23T04:09:19Z

Just a doubt. I add a header in t.txt which do not appear in new.txt.

PaulWessel · 2022-01-23T04:12:42Z

OK, will have a look.

PaulWessel · 2022-01-23T04:23:33Z

Hm, I added both a header (# My header) and a segment header (> My segment):

gmt spatial t.txt -W60n/100k 
> My segment
-1.33158486117	0.945529552785
-0.5	1.5
0.25	2
0.75	2
1.25	3

Usually headers are not preserved on output. Adding -h gives

# My Header
# Command : gmtspatial t.txt -W60n/100k -h
> My segment
-1.33158486117	0.945529552785
-0.5	1.5
0.25	2
0.75	2
1.25	3
1.65260912694	3.804522831

What is it that is not working for you?

Esteban82 · 2022-01-23T04:27:51Z

Ok, sorry. Perfect.

joa-quim · 2022-01-23T10:57:10Z

src/gmtspatial.c

 	}
 #endif

+	if (Ctrl->W.active) {	/* Calculate centroid and polygon areas or line lengths and place in segment headers */


Is this comment correct?

joa-quim · 2022-01-23T11:05:01Z

Would be better to slow down a bit the merging step. A couple of things.

The docs should mention if this calculation is spherical or ellipsoidal
It would preferable that the tests were numeric instead of graphical
Allow also to specify the azimuth along which de extrapolation is done. That azimuth could even be calculated from a best fit of all points in the segment. Useful if there is some noise in data.

PaulWessel · 2022-01-23T19:57:08Z

Those are good points, @joa-quim . I think @Esteban82 was eager to try this for his work. I guess a good general guideline is to let the one who opened the PR get to close it since that person may discover minor stuff that should be added to the PR before merging.

PaulWessel · 2022-01-23T20:18:23Z

I have fixed the dumb copy/paste comment and added explanation to the docs that the geographic calculations are spherical. I could code up Karney 2012 but this is already in proj and/or GDAL so we could add this as needed. It makes no difference for short extensions but will show up for long extensions.

As for the azimuthal options. If these are sensible and useful they can be added. I did think of azimuths but best fit on a sphere is what fitcircle does for instance, but it does not report azimuths so we would have to add some coding here. And again this is pseudo-elliptical (authalic parameters for a sphere). For discussion, I can imagine the new modifier

+a[azim[/azim2]]

where +a (no args) means determine azimuth from best fit great/small circle (or least-squares line for Cartesian), +aazim means use that azimuth at either end (of the ends selected if +f|l is used) and +aazimf/asiml for separate azimuths. Note that for long lines and lines close to poles the single azimuth option would give very different results at either end.

Might it be useful to have modifiers for adding a deviation to the selected azimuth? E.g., +a+d10/5 would add 10 degrees to the azimuth at start and 5 at the end?

joa-quim · 2022-01-23T20:30:18Z

Might it be useful to have modifiers for adding a deviation to the selected azimuth? E.g., +a+d10/5 would add 10 degrees to the azimuth at start and 5 at the end?

Seems adding too much complexity to syntax for very little optional gain.

Karney 2012 but this is already in proj and/or GDAL

Yes, it's in proj and that's what I used in my Julia example. I think using proj would have been better as it gives ellipsoidal instead of spherical. And about this I made a little test this morning that shows that you are using the azimuth at the start instead of end point. I'm 100% clear on this but I think we should use the azimuth at the end point.

 julia> d, az1, az2 = invgeod([-66.04971205 -55.14285297], [-65.90790703 -55.03662696], proj="+proj=longlat +ellps=sphere")
([14864.354381839295], [37.43707316842508], [37.3207859763504])

# Use azim at first point
julia> geod([-65.90790703 -55.03662696], 37.43707316842508, 15000, proj="+proj=longlat +ellps=sphere")
([-65.76518956299792 -54.92943115095463], 37.32019004431618)

# Use azim at second point
julia> geod([-65.90790703 -55.03662696], 37.3207859763504, 15000, proj="+proj=longlat +ellps=sphere")
([-65.7655687941153 -54.929265381988785], 37.204213554034325)

and the GMT solution

 gmtspatial estica.dat -W15000e+l
-66.04971205    -55.14285297
-65.90790703    -55.03662696
-65.7651897905  -54.9294313224

joa-quim · 2022-01-23T20:38:54Z

Yet another useful option would be +nndists to turn the dist into an increment and thus compute points along the line ... which at the end corresponds to compute a chunk of a great circle, so perhaps there are better ways to do that.
Ghrr, it's not a great circle. These things get confusing very rapidly.

PaulWessel · 2022-01-23T20:43:02Z

Perhaps let that happen in project -G instead. I will check the azimuths, the intent was to get the azimuth at the end points pointing outwards.

PaulWessel · 2022-01-24T01:26:22Z

Here is the result for geodesic (Vicenty):

gmt spatial estica.dat-W15k -je --FORMAT_FLOAT_OUT=%.16g
-66.19291065258025	-55.24983297958808
-66.04971205	-55.14285297
-65.90790703	-55.03662696
-65.76547098032331	-54.92964199505121

So for the end point this differs from your Kerney solution by 1.6e-13 in longitude and 2e-14 in latitude, i.e., ~9 nm. However, at the first point I get 25 meters difference, so we need to resolve this. I am doing this:

To extend the start point backwards I compute the back-azimuth from the 2nd point to the 1st point; this is then the azimuth of extending the line further: par[0] = gmt_az_backaz (GMT, S->data[GMT_X][1], S->data[GMT_Y][1], S->data[GMT_X][2], S->data[GMT_Y][2], true);
To extend the end point forward I compute the azimuth from the 1st point to the 2nd point; this is then the azimuth of extending the line at the end: par[0] = gmt_az_backaz (GMT, S->data[GMT_X][S->n_rows-3], S->data[GMT_Y][S->n_rows-3], S->data[GMT_X][S->n_rows-2], S->data[GMT_Y][S->n_rows-2], false);

Code says:

	**	latS, lonS -- latitude and longitude of first point in radians.
	**	latE, lonE -- latitude and longitude of second point in radians.
	**
	**	OUTPUT
	**  	faz -- azimuth from first point to second in radians clockwise from North.
	**	baz -- azimuth from second point back to first point.
	**	bool back_az controls which is returned

I think this is right, no? Note that my row indices above is after I have added the blank rows and start and end, so the two input points are at index 1 and 2, and also n-3 and n-2 (since n = 4).

joa-quim · 2022-01-24T01:41:13Z

What I do is to use baz+180 and for the first point I swapped 1->2 to 2->1

So it seems the same procedure. But 20-40 differences is what I found when comparing the effect of using faz and baz+180 (tried it only at end point).

Note, you can use the geod command to get the same as I got with the Julia wrapped functions. It's a bit annoying as it wants lat lon instead of the civilized lon lat
Using faz in this example

geod +ellps=WGS84 -I -f %.15f -
-55.14285297 -66.04971205 -55.03662696 -65.90790703
37.498033893416832      -142.618253298774135    14893.889

geod +ellps=WGS84 -f %.15f -
-55.03662696 -65.90790703 37.498033893416832
-55.036626960000007     -65.907907030000004     -142.501966106583154

PaulWessel · 2022-01-24T02:03:09Z

Note that usually the azimuths at one end + 180 is not the same as at the other end:

gmt mapproject -AB -je t.txt -fg
-66.04971205	-55.14285297	NaN
-65.90790703	-55.03662696	-142.618253299
gmt mapproject -AF -je t.txt -fg
-66.04971205	-55.14285297	NaN
-65.90790703	-55.03662696	37.4980338934
gmt math -Q 37.4980338934 180 SUB =
-142.501966107

-142.501966107 ~= -142.618253299

PaulWessel · 2022-01-24T02:15:13Z

Note: As the two points converge this difference goes to zero of course, but since the line segment is finite and a small section of the great circle (or geodesic), the azimuth is a continuous function along the geodesic and cannot be the same at two different points.

joa-quim · 2022-01-24T02:19:19Z

I know. I'm adding 180 when it gives the baz to get a forward azim at point 2 like the example I showed when using the geod exe.

PaulWessel · 2022-01-24T02:19:22Z

So I think my implementation is correct. The baz = az+180 only applies at the same point, not at the end points. Try with the line segment given by (0,0) and (45,45) and the difference in the two azimuths is huge: 35.4100589051 vs -125.109226172 and if you subtract 180 from 35.4100589051 you get -144.589941095.

joa-quim · 2022-01-24T02:22:16Z

In fact I'm not adding any 180 when using Julia because what the geod function returns is the forwards azims at point 1 and point 2. Is the geod cli that returns a backward azim at point 2. Then if one wants to extend from there (pt2) we must add 180.

joa-quim · 2022-01-24T02:29:42Z

OK, take that example. What I do is
54.890773827875094 = -125.109226172124906 + 180
where -125.109226172124906 is the baz at 45,45

geod +ellps=WGS84 -I -f %.15f -
0 0 45 45
35.410058905065931      -125.109226172124906    6662472.718

geod +ellps=WGS84 -f %.15f -
45 45 54.890773827875094 10000
45.051705283658144      45.103846144750108      -125.035762716108209

PaulWessel · 2022-01-24T02:31:13Z

Please do the case for the line (0,60) - (30,85) [two points with 100k extension at either ends so it is easier to check for crazy results. I can then compare

PaulWessel · 2022-01-24T02:37:23Z

You are right - I need to flip the calls in gmtspatial. The az_backaz function is OK but I was using it wrong and the extend.sh test was too tiny distances to tell.

joa-quim · 2022-01-24T02:43:33Z

Ok but since I did (in Julia because the geod cli (you also have it) is very boring)

julia> d, az1, az2 = invgeod([0. 60], [30 85])
([2.876587608896881e6], [5.756074039070723], [35.09264842252227])

julia> geod([30 85], 35.09264842252227, 100000)
([36.88394599688237 85.70175392814406], 41.954127083279346)

julia> d, az1, az2 = invgeod([30. 85], [0. 60])
([2.876587608896881e6], [-144.90735157747773], [-174.2439259609293])

julia> geod([0. 60], -174.2439259609293, 100000)
([-0.17503387956791308 59.1067786323064], -174.39482758940514)

Let gmtspatial -W extend segment

e04d513

This PR follows the comments on the forum post and adds a new -W option that is used to extend existing lines in one or both directions by a given amount. A test has been added and the documentation has been updated.

PaulWessel added the new core module feature PR that implements a new core module feature label Jan 23, 2022

PaulWessel added this to the 6.4.0 milestone Jan 23, 2022

PaulWessel requested review from Esteban82, joa-quim and maxrjones January 23, 2022 03:15

PaulWessel self-assigned this Jan 23, 2022

Esteban82 approved these changes Jan 23, 2022

View reviewed changes

Esteban82 merged commit 12843ce into master Jan 23, 2022

Esteban82 deleted the extend-segments branch January 23, 2022 03:53

joa-quim reviewed Jan 23, 2022

View reviewed changes

maxrjones added the add-changelog Add PR to the changelog label Jan 25, 2022

Let gmtspatial -W extend segments from one of both ends #6255

Let gmtspatial -W extend segments from one of both ends #6255

Uh oh!

Conversation

PaulWessel commented Jan 23, 2022

Uh oh!

Esteban82 left a comment

Choose a reason for hiding this comment

Uh oh!

Esteban82 commented Jan 23, 2022

Uh oh!

PaulWessel commented Jan 23, 2022

Uh oh!

PaulWessel commented Jan 23, 2022

Uh oh!

Esteban82 commented Jan 23, 2022

Uh oh!

joa-quim Jan 23, 2022

Choose a reason for hiding this comment

Uh oh!

PaulWessel Jan 23, 2022

Choose a reason for hiding this comment

Uh oh!

joa-quim commented Jan 23, 2022

Uh oh!

PaulWessel commented Jan 23, 2022

Uh oh!

PaulWessel commented Jan 23, 2022

Uh oh!

joa-quim commented Jan 23, 2022

Uh oh!

joa-quim commented Jan 23, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

PaulWessel commented Jan 23, 2022

Uh oh!

PaulWessel commented Jan 24, 2022

Uh oh!

joa-quim commented Jan 24, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

PaulWessel commented Jan 24, 2022

Uh oh!

PaulWessel commented Jan 24, 2022

Uh oh!

joa-quim commented Jan 24, 2022

Uh oh!

PaulWessel commented Jan 24, 2022

Uh oh!

joa-quim commented Jan 24, 2022

Uh oh!

joa-quim commented Jan 24, 2022

Uh oh!

PaulWessel commented Jan 24, 2022

Uh oh!

PaulWessel commented Jan 24, 2022

Uh oh!

joa-quim commented Jan 24, 2022

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

joa-quim commented Jan 23, 2022 •

edited

Loading

joa-quim commented Jan 24, 2022 •

edited

Loading