How to use gmt spatial -Sb correctly?

I need to calculate some buffer polygons around some lines. According to the docs, gmt spatial -Sb should to the trick:

-Sbwidth which computes a buffer polygon around lines

Source

However it did not exactly do what I expected. I made three test cases to figure out what’s going on:

Green is the original line, red is the computed buffer polygon. Test script is at the end of this post. Maybe I’m making some large mistakes but I observed the following problems:

  • General:
    • The scaling of the units seems to be way off. I wanted a 10 nautical mile (NM) buffer around the line but this led to a polygon way outside the plot area. had to use 0.05 NM to get the present result.
    • Different scaling in N-S and E-W directions. The buffer seems smaller in E-W direction compared to the buffer in N-S direction. Maybe a bug where it wasn’t taken into account that 1’ latitude isn’t 1 NM except on the equator? This Lambert conic conformal projection should make the buffer almost identical in both directions.
  • Test 1 (upper left) did almost what I expected, however there are two stray lines pointing to the west?
  • Test 2 (middle) as expected except for the points written under General.
  • Test 3 (lower right) has similar problems as Test 1 (stray lines) and additionally the buffer simply bridges the gap and doesn’t follow the green line on the inside. I expected something like Test 2 but with the bends.

As said, I’m not sure if I did everything the proper way. Your ideas and help is greatly appreciated. Am I making some major mistakes somewhere?

Using GMT version: 6.2.0_43ff8ef_2021.05.25

Test script to play with:

#!/usr/bin/env bash

echo "Using GMT version: $(gmt --version)"

cat > testline1.txt << EOF
>
7.8 50.4
7.8 50.6
8.1 50.6
8.1 50.4
7.8 50.4
EOF

cat > testline2.txt << EOF
>
8   50
9.5 50.5
EOF

cat > testline3.txt << EOF
>
9   49.6
9   50
9.8 50
9.8 49.8
9.4 49.8
9.4 49.6
EOF

gmt spatial testline1.txt -Sb0.05n > testline1_buffer.txt
gmt spatial testline2.txt -Sb0.05n > testline2_buffer.txt
gmt spatial testline3.txt -Sb0.05n > testline3_buffer.txt

gmt begin buffer png
  gmt basemap -R007:30/010:00/49:30/50:50 \
    -JL008:34:13.64/50:01:59.90/49:57/50:23/15c -B
    
  gmt plot testline1.txt -W1p,green
  gmt plot testline1_buffer.txt -W1p,red

  gmt plot testline2.txt -W1p,green
  gmt plot testline2_buffer.txt -W1p,red

  gmt plot testline3.txt -W1p,green
  gmt plot testline3_buffer.txt -W1p,red
gmt end show

The stray lines are clearly a bug. The rest is tougher. I recommend reading a bit of this thread. And, BTW, I think Mirone does it nicely.

Precisely

Good evening @Joaquim, I assume you did the plot in Mirone? What kind of projection did you use? Mercator?

Yes, it’s a Mirone plot but the projection if it has a name I don’t know it. It’s a linear one that has a 1:1 scale at middle latitude. But that is not the real point. If you have read the mails I linked you saw that the main problem is that GEOS lib (the guy who computes the buffers) works only for Cartesians and the mail discuss strategies when using lon,lat data (basically the advice is to project to an equal-area, compute the buffer and project back to geogs. But this is not implemented in gmtspatial.

In Mirone I programmed another algorithm, (based on the Matlab one) where each vertex is the source of a geographical circle (radius = buffer width) and the buffer is the union of all circles and rectangles connecting not touching circles. This is not completely right as I should interpolate the straight lines with a 2*buff_width such that vertex insured circles touch each other. But for my usage at the time, that was good enough.

In Mercator

Yes I read the mailing list thread. Thank you for the link. I’m not sure if the equal area work-around is good enough for higher latitudes. I hoped with this experiment to find another approach for my problems from this thread or at least to get a nice outline for visualisation purposes.

Perhaps make this into a bug report on GitHub? Also, if -Sb requires Cartesian data then this should be pointed out in the docs.

Issue opened:

@KristofKoch With new buffergeo function in Julia (master) you can now

Dk = buffergeo("tobuff.dat", width=5, unit=:Nautical);
imshow(Dk)

Neat! I’ll give it a test drive right away. But also limited to Cartesian data?

On the contrary, this one is purely geographic (spherical or ellipsoidal [default]) and believing in this Esri site not even ArcGIS has one (theirs is only for points or multi-points), whilst QGis uses the goto-equal-area-n-return and R does it only with (visible) rectangles.

Matlab also has one that is spherical only.

In this fig the red line represents the Matlab solution. I’m still not clear why they differ so much on the polar region. I interpolate the line segments along great circles, Matlab rotates the segments to put them on the equator and then computes the buffer along small circles of latitude = to buffer width so I would say that they also interpolate the line segment along great circles but it looks more as if the interpolation is done along rumb lines.

Wow, that’s impressive @Joaquim!

When trying to recreate it I only get error messages. I installed julia via github and have now version 1.7.0-DEV.1268 (2021-06-07). For troubleshooting I went back to copy & pasting examples from the GMT.jl docs but no luck either:

julia> coast(region=:g, proj=(name=:ortho, center=(-75,41)), frame=:g, res=:crude, area=5000,
             land=:pink, ocean=:thistle, figsize=10, show=true)
ERROR: UndefVarError: coast not defined
Stacktrace:
 [1] top-level scope
   @ REPL[9]:1

julia> 

Do you know what I’m doing wrong here? I assume it is some kind of beginner mistake.

You are forgetting to do

using GMT

before the commands.

Hi @Joaquim this was definitly one of my errors. However I still get errors:

julia worldsnake.jl                                                                                   
ERROR: LoadError: could not load symbol "pj_get_spheroid_defn":
/usr/lib/libproj.so.22: undefined symbol: pj_get_spheroid_defn
Stacktrace:
 [1] pj_get_spheroid_defn
   @ ~/.julia/packages/GMT/Ixkua/src/proj_utils.jl:416 [inlined]
 [2] get_ellipsoid(projPJ_ptr::Ptr{Nothing})
   @ GMT ~/.julia/packages/GMT/Ixkua/src/proj_utils.jl:390
 [3] invgeod(lonlat1::Matrix{Float64}, lonlat2::Matrix{Float64}; proj::String, s_srs::String, epsg::Int64)
   @ GMT ~/.julia/packages/GMT/Ixkua/src/proj_utils.jl:126
 [4] buffergeo(line::Matrix{Float64}; width::Int64, unit::Symbol, np::Int64, flatstart::Bool, flatend::Bool, proj::String, epsg::Int64, tol::Float64)
   @ GMT ~/.julia/packages/GMT/Ixkua/src/proj_utils.jl:241
 [5] buffergeo(D::Vector{GMT.GMTdataset}; width::Int64, unit::Symbol, np::Int64, flatstart::Bool, flatend::Bool, epsg::Int64, tol::Float64)
   @ GMT ~/.julia/packages/GMT/Ixkua/src/proj_utils.jl:232
 [6] buffergeo(fname::String; width::Int64, unit::Symbol, np::Int64, flatstart::Bool, flatend::Bool, epsg::Int64, tol::Float64)
   @ GMT ~/.julia/packages/GMT/Ixkua/src/proj_utils.jl:225
 [7] top-level scope
   @ ~/Dokumente/GMTjl/worldsnake.jl:6
in expression starting at /home/kristof/Dokumente/GMTjl/worldsnake.jl:6

Would you mind you kindly sharing your example code? Maybe I can figure it out myself.

Hmmm, that’s worse. It means that your proj4 library is old or for some reason it misses the pj_get_spheroid_defn function.
At some point in the development I saw that error too (don’t remember anymore where)
It works in the 3 Linux’s that I have access.

A very simple command to plot a buffer is

imshow(buffergeo([0 0; 5 5], width=100, unit=:k))

What does this print?

GMT.proj_info()

I tried to update my julia installation via git pull followed by make clean && make as recommended in the docs but now the make process fails:

[…]
    CC src/llvm-demote-float16.o
    LINK usr/lib/libjulia-internal.so.1.7
/usr/bin/ld: ./codegen.o: in function `__gnu_cxx::new_allocator<llvm::Type*>::allocate(unsigned long, void const*)':
/usr/include/c++/11.1.0/ext/new_allocator.h:110: undefined reference to `std::__throw_bad_array_new_length()'
/usr/bin/ld: ./llvm-late-gc-lowering.o: in function `__gnu_cxx::new_allocator<int>::allocate(unsigned long, void const*)':
/usr/include/c++/11.1.0/ext/new_allocator.h:110: undefined reference to `std::__throw_bad_array_new_length()'
/usr/bin/ld: ./llvm-late-gc-lowering.o: in function `__gnu_cxx::new_allocator<unsigned int>::allocate(unsigned long, void const*)':
/usr/include/c++/11.1.0/ext/new_allocator.h:110: undefined reference to `std::__throw_bad_array_new_length()'
/usr/bin/ld: ./llvm-late-gc-lowering.o: in function `__gnu_cxx::new_allocator<int>::allocate(unsigned long, void const*)':
/usr/include/c++/11.1.0/ext/new_allocator.h:110: undefined reference to `std::__throw_bad_array_new_length()'
/usr/bin/ld: /usr/include/c++/11.1.0/ext/new_allocator.h:110: undefined reference to `std::__throw_bad_array_new_length()'
/usr/bin/ld: ./llvm-late-gc-lowering.o:/usr/include/c++/11.1.0/ext/new_allocator.h:110: more undefined references to `std::__throw_bad_array_new_length()' follow
collect2: Fehler: ld gab 1 als Ende-Status zurück
make[1]: *** [Makefile:300: /home/kristof/Dokumente/github/julia/usr/lib/libjulia-internal.so.1.7] Fehler 1
make: *** [Makefile:76: julia-src-release] Fehler 2

I feel I wrecked something there …

The problem is not Julia. Just use the binaries they provide. 1.6 or 1.7-DEV are OK. The problem is GDAL+PROJ4 (and linux). Again, what does it says?

Ok, I reverted back to julia 1.6.1. Here you go:

julia> GMT.proj_info()
Rel. 8.0.1, March 5th, 2021
Located at: /home/kristof/.local/share/proj:/usr/share/proj:/usr/share/proj