Midpoint of a line

EvanPi · July 1, 2025, 8:23pm

Hi,
I am trying to get a midpoint of a line.
awk -F, '{mid_lon=($1+$3)/2; mid_lat=($2+$4)/2; print mid_lon, mid_lat}' dredgeinfo.txt > dredgemidpoint.txt
that should take the average of 2 points, but it appears that this point isnt quite on line conencting the two points. i have a feeling this is some sort of projection thing, where the midpoint isnt calculated along a great circle or something similar. How can I do this with GMT?

Location of 2 endpoints:
-104.888 -4.559
-104.877 -4.553

thwbecker · July 1, 2025, 9:44pm

GMT by default uses a great circles to connect geographic points unless that is suppressed with the -A flag. Turned around, to get the mid point distance in gegraophic coordinates you can use a script like this one

github.com/thwbecker/awk

pointfromtwoloc.awk

main

#
# given two lon lat locations, compute point at relative distance along 
# great circle
#
# input:
# lon1 lat1 lon2 lat2 frac [other_stuff]
#

BEGIN{
    pi = 3.141592653589793238;
    twopi = 2.0*pi;
    f = pi/180;
}
{
    if((substr($1,1,1)!="#")&&(NF>=5)){
	lon1 = $1;
	lat1 = $2;
	lon2 = $3;
	lat2 = $4;
	frac = $5;

This file has been truncated. show original

assuming spherical Earth or use a projection package to get more accurate results.

jaltekruse · July 2, 2025, 5:17pm

Your awk print statement may be truncating the coordinates, try printf("%f %f",mid_lon,mid_lat) (when I tried this, I get the same results from the GMT solution below and both plot on the line). To stay within GMT, I think gmt spatial can find the midpoint:

$ cat tmp_line.gmt 
>
-104.888 -4.559
-104.877 -4.553

$ gmt spatial -Qk tmp_line.gmt
-104.8825       -4.556  1.38849721954

Here, gmt spatial prints the coordinates of the midpoint and the total length of the line (1.39 km).

mapproject can calculate the distance and azimuths from the endpoints to the midpoint to verify that the azimuths from the endpoints to the midpoint are 180 degrees apart and the distances from the endpoints to the midpoint are half of the total length:

$ gmt mapproject -G-104.8825/-4.556+uk tmp_line.gmt -Af-104.8825/-4.556
>
-104.888        -4.559  241.312977761   0.694247508686
-104.877        -4.553  61.3135153114   0.694249710865

Joaquim · July 2, 2025, 10:05pm

Well, strictly speaking it doesn’t compute the mid-point. It does compute the centroid of the polygon. But since the it closes the line and apparently thinks it’s a polygon, the result is the same

mkononets · July 4, 2025, 10:09pm

jaltekruse:

$ gmt mapproject -G-104.8825/-4.556+uk tmp_line.gmt -Af-104.8825/-4.556
>
-104.888        -4.559  241.312977761   0.694247508686
-104.877        -4.553  61.3135153114   0.694249710865

also @Joaquim

Something I got VERY puzzled over: when I run the above without -G..., the calculate azimuths become different:

gmt mapproject  tmp_line.gmt -Af-104.8825/-4.556
-104.888	-4.559	241.389540334
-104.877	-4.553	61.389540334

i.e. 61.389540334 vs 61.3135153114

Is this me doing something wrong? Thanks! I’ve been making some azimuth calculations for oblique Mercator -Joa... specification, and this makes a lot of difference.

Joaquim · July 5, 2025, 10:49am

OK, things I found.
First let me start with the correct ellipsoidal answer to the mid point (these GMT.jl functions use PROJ).

julia> invgeod([-104.888, -4.559], [-104.877, -4.553])
(1389.3356054500066, 61.47443481564943, 61.473561047920725)

we get, distance, azimuth at starting point, azimuth at end point. Now compute the midpoint

julia> geod([-104.888, -4.559], 61.47443481564943, 1389.3356054500066/2)
([-104.88249997720496 -4.556000021167548], 61.47399778643971)

second value is the azimuth at midpoint. The midpoint coordinates are a bit different of those by gmtspatial. Likely the reason is that gmtspatial is not able to compute centroids using the ellipsoidal coordinates.

Now mapproject -A has the same problem. It ignores the ellipsoidal request and the azimuths are different then those of GMT.jl

C:\v>mapproject a.dat  -Af-104.88249997720496/-4.556000021167548 -je
-104.888        -4.559  241.389810108
-104.877        -4.553  61.3892705594

C:\v>mapproject a.dat  -Af-104.88249997720496/-4.556000021167548
-104.888        -4.559  241.389810108
-104.877        -4.553  61.3892705594

This is a bug, Specially since if we use the -G and -je then results match those of PROJ … and we get a weird warning

mapproject a.dat -G-104.88249997720496/-4.556000021167548 -Ab-104.88249997720496/-4.556000021167548 -je
mapproject [WARNING]: Your distance mode (Great Circle) differs from your -j option (Geodesic) which takes precedence.
mapproject [WARNING]: Your distance mode (Great Circle) differs from your -j option (Geodesic) which takes precedence.
-104.888        -4.559  61.4744348127   694.667802659
-104.877        -4.553  -118.526438955  694.667802659

Now both GMT and PROJ azimuths (61.4744348127 vs 61.47443481564943) match up to the 8th decimal. Note that I used -Ab instead of -Af given the bit weird way that GMT uses to define with respect to which point we refer as forward and backward azimuths.

Joaquim · July 5, 2025, 12:00pm

More. The problem is that with only -A it is not detecting that the input is geographical. Adding -fg makes it work.

mapproject a.dat  -AB-104.88249997720496/-4.556000021167548 -fg
-104.888        -4.559  61.4744348127
-104.877        -4.553  -118.526438955

and

mapproject a.dat  -Ab-104.88249997720496/-4.556000021167548 -fg -je
mapproject [WARNING]: Your distance mode (Great Circle) differs from your -j option (Geodesic) which takes precedence.
-104.888        -4.559  61.4744348127
-104.877        -4.553  -118.526438955

These upper case F,B must be some old stuff before the -je introduction and should be removed in favor of later.

mkononets · July 5, 2025, 12:24pm

Thanks a lot!

what about -je and gmt spatial?

gmt spatial -Qk tmp_line.gmt -je
gmtspatial [WARNING]: Your distance mode (Great Circle) differs from your -j option (Geodesic) which takes precedence.
-104.8825	-4.556	1.38933560532
gmt spatial -Qk tmp_line.gmt
-104.8825	-4.556	1.38849721954

mkononets · July 5, 2025, 12:35pm

that’s kind of weird as CLI geod reports -118.52643895207179 instead of 61.473561047920725

echo -4.559 -104.888 -4.553 -104.877 | invgeod +ellps=WGS84 -F "%.17g" -f "%.17g" 
61.474434815656934	-118.52643895207179	1389.3356054500071

Joaquim · July 5, 2025, 1:49pm

Must see what invgeod CLI says it reports. The one via GMT.jl says

Returns

dist - A scalar with the distance between point 1 and point 2 (meters). Or a vector when lonlat1|2 have more than one pair of points.

az1 - azimuth at point 1 (degrees) ∈ [-180, 180)

az2 - (forward) azimuth at point 2 (degrees) ∈ [-180, 180)

Joaquim · July 5, 2025, 1:50pm

That’s what I said. Centroid computation seem to ignore -j (but length, not).

mkononets · July 5, 2025, 2:53pm

CLI geod apparently reports the same, forward and back azimuths and distance between
an initial and terminus point latitudes and longitudes (inverse). 180 - 118.52643895207179 = 61.47356104792821 is the same as 61.473561047920725 to the 11th decimal digit.

Both man geod and geod — PROJ 9.6.2 documentation are less detailed though:

geod (direct) and invgeod (inverse)  perform  geodesic  (Great  Circle)
       computations  for determining latitude, longitude and back azimuth of a
       terminus point given a initial point latitude, longitude,  azimuth  and
       distance (direct) or the forward and back azimuths and distance between
       an initial and terminus point latitudes and longitudes (inverse).   The
       results  are  accurate to round off for |f| < 1/50, where f is flatten‐
       ing.
...
       -p     This option causes the azimuthal values to be output as unsigned
              DMS numbers between 0 and 360 degrees. Also note -f.

Thanks a lot anyway for reminding that geod exists and does its thing.

Joaquim · July 5, 2025, 2:58pm

Not exactly. It says and back azimuths, whilst on GMT.jl they are both forward, but backwards are also possible if backward=true is set.

dist, az1, az2 = invgeod(lonlat1::Vector{<:Real}, lonlat2::Vector{<:Real}; proj::String="",
                         s_srs::String="", epsg::Integer=0, backward=false)

Joaquim · July 5, 2025, 3:17pm

Fixes submitted in mapproject -A does not respect ellipsoidal request · Issue #8759 · GenericMappingTools/gmt · GitHub

mkononets · July 5, 2025, 3:36pm

one more question about gmt spatial vs geod:

Joaquim:

julia> geod([-104.888, -4.559], 61.47443481564943, 1389.3356054500066/2)
([-104.88249997720496 -4.556000021167548], 61.47399778643971)
second value is the azimuth at midpoint. The midpoint coordinates are a bit different of those by gmtspatial. Likely the reason is that gmtspatial is not able to compute centroids using the ellipsoidal coordinates.

gmt spatial -Qk tmp_line.gmt -je --FORMAT_FLOAT_OUT=%.17g
gmtspatial [WARNING]: Your distance mode (Great Circle) differs from your -j option (Geodesic) which takes precedence.
-104.88249999999999	-4.556	1.3893356053178714

Actually, for the centroid coordinates the difference is very small, see below.

longitude: -104.88249999999999-(-104.88249997720496) = -0.00000002279503
latitude: -4.556-(-4.556000021167548) = 0.000000021167548
that’s ~ 2.3e-8, the 8th decimal.
With the distance 1.3893356053178714 - 1389.3356054500066/1000 = -0.0000000001321352
that’s 10th decimal digit if the distance is in km or 7th if in meters.

Is the above not negligible? Would be good to know sort of boundary when one needs to start investigating the discrepancies.

Joaquim · July 5, 2025, 3:47pm

Confess that I didn’t pay much attention to the gmtspatial case. Yes, those differences are negligible. For the centroid/midpoint I would expect that they are due to different methods used to compute it, but for the length I feel some tickles on why they don’t agree to more decimals as they are using (I think) the same equations, but nothing that I want to time in investigating.

mkononets · July 5, 2025, 3:52pm

Thanks anyway!

mkononets · July 9, 2025, 8:13pm

@Joaquim
Sorry for bothering with this again, but I think mapproject's -W option does not respect -je either.

I am not sure it is supposed to in this situation, but asking usually does not harm.

I tried to use mapproject -Joa.../-Job... -We... to calculate a rectangle of an oblique area, like it is shown on the figure in mapproject’s documentation, left panel.

Here I am using the line from the topic, the ends and the calculated midpoint. The midpoint is needed for -Joa/-Job... specification. Then a rectangle 100 m (0.1 km) wide is specified as the oblique rectangular region around the line. Half-length is also needed for the correct -R... specification.

For this I am calling mapproject like below:

gmt mapproject -Job-104.88249997720496/-4.5560000211675495/-104.877/-4.553/3i -R0.69424861297385021/0.100+uk -We+n2 -je
255.112438522	-4.5597859581
255.1234276	-4.55379191313
255.12256147	-4.55221404886
255.111572409	-4.55820808067
255.112438522	-4.5597859581

gmt mapproject -Job-104.88249997720496/-4.5560000211675495/-104.877/-4.553/3i -R0.69424861297385021/0.100+uk -We+n2
255.112438522	-4.5597859581
255.1234276	-4.55379191313
255.12256147	-4.55221404886
255.111572409	-4.55820808067
255.112438522	-4.5597859581

^^^ output, i.e. rectangle coords with and without -je match. The same is true for the -Joa... case (not shown). For mapproject -Joa... I calculated azimuth externally, see more below.

My posting should have probably ended up here, but I could not just stop writing.

Another part of the same problem (I guess) was when I compared azimuth-based equator -Joa... -We+n2 (using separate azimuth calculations mapproject -A... -G...) vs two-point-based equator -Job... -We+n2

When I used spherical azimuth from mapproject -G... -A... (without -je), the rectangles calculated using mapproject -Joa... -We+n2 and mapproject -Job... -We+n2, matched. Please see Spaghetti 1 listing below.

When I used elliptical azimuth from mapproject -G... -A..., the resulting rectangles were visibly different. Please see Spaghetti 2 listing below.

Spaghetti 1, spherical mapproject azimuth, the two oblique rectangles match well:

geod +ellps=WGS84 -F %.17g -f %.17g +lat_1=-4.559 +lon_1=-104.888 +lat_2=-4.553 +lon_2=-104.877 +n_S=2
-4.5590000000000002	-104.88800000000001
-4.5560000211675495	-104.88249997720496
-4.5529999999999999	-104.877

gmt mapproject -G-104.88249997720496/-4.5560000211675495+uk -Af-104.88249997720496/-4.5560000211675495 --FORMAT_FLOAT_OUT=%.17g
-104.877	-4.5529999999999999	61.313245074112395	0.69424861297385021

#gmt mapproject Azimuth is 61.313245074112395, LENGTH_KM is 0.69424861297385021

gmt mapproject -Job-104.88249997720496/-4.5560000211675495/-104.877/-4.553/3i -R0.69424861297385021/0.100+uk -We+n2
255.112438522	-4.5597859581
255.1234276	-4.55379191313
255.12256147	-4.55221404886
255.111572409	-4.55820808067
255.112438522	-4.5597859581

gmt mapproject -Joa-104.88249997720496/-4.5560000211675495/61.313245074112395/3i -R0.69424861297385021/0.100+uk -We+n2
255.112438522	-4.5597859581
255.1234276	-4.55379191307
255.12256147	-4.5522140488
255.111572409	-4.55820808067
255.112438522	-4.5597859581

^^^ the two oblique rectangles match well

Spaghetti 2, the elliptical azimuth using mapproject -G... -A... -je, the two oblique rectangles at the end do not match:

geod +ellps=WGS84 -F %.17g -f %.17g +lat_1=-4.559 +lon_1=-104.888 +lat_2=-4.553 +lon_2=-104.877 +n_S=2
-4.5590000000000002	-104.88800000000001
-4.5560000211675495	-104.88249997720496
-4.5529999999999999	-104.877

#LON_CENTER is -104.88249997720496, LAT_CENTER is -4.5560000211675495

gmt mapproject -G-104.88249997720496/-4.5560000211675495+uk -Af-104.88249997720496/-4.5560000211675495 --FORMAT_FLOAT_OUT=%.17g -je
mapproject [WARNING]: Your distance mode (Great Circle) differs from your -j option (Geodesic) which takes precedence.
mapproject [WARNING]: Your distance mode (Great Circle) differs from your -j option (Geodesic) which takes precedence.
-104.877	-4.5529999999999999	61.473997783466466	0.69466780265863348

#gmt mapproject Azimuth is 61.473997783466466, LENGTH_KM is 0.69466780265863348

# geod azimuth calculations for comparison
echo -4.5560000211675495 -104.88249997720496 -4.553 -104.877 | geod +ellps=WGS84 -F %.17g -f %.17g -I
61.473997786429415	-118.52643895208958	694.66780272511619

#geod Azimuth is 61.473997786429415, LENGTH_M is 694.66780272511619

gmt mapproject -Job-104.88249997720496/-4.5560000211675495/-104.877/-4.553/3i -R0.69466780265863348/0.100+uk -We+n2
255.112435205	-4.55978776767
255.123430918	-4.55379010348
255.122564788	-4.55221223922
255.111569091	-4.55820989024
255.112435205	-4.55978776767

gmt mapproject -Joa-104.88249997720496/-4.5560000211675495/61.473997783466466/3i -R0.69466780265863348/0.100+uk -We+n2
255.112424564	-4.5597735875
255.123437115	-4.55380669965
255.122575429	-4.55222641921
255.111562895	-4.55819329394
255.112424564	-4.5597735875

^^^ the two oblique rectangles at the end do not match

Joaquim · July 9, 2025, 11:57pm

A very quick comment without having read your questions with the needed attention.
I don’t think that Oblique Marcator projection have an ellipsoidal version. Only spherical. This would explain at least part of your equalities.

mkononets · July 10, 2025, 1:33pm

That would explain everything I found, the fact that Oblique Mercator projection is only available as spherical. That is, only the spherical distances and azimuths produce the matching oblique rectangles using mapproject ...-We+n2.

Thanks again!