Feature #1136

grdcontour -Q: specify length instead of number of points

Added by Andreas about 4 years ago. Updated over 3 years ago.

Status:ClosedStart date:2017-07-26
Priority:NormalDue date:
Assignee:Paul% Done:


Target version:Candidate for next minor release


It looks like grdcontour does a contour length calculation and compares it to a minimum value under the hood (from grdcontour.c);

#define GRDCONTOUR_MIN_LENGTH 0.01        /* Contours shorter than this are skipped */

Contours shorter than a minimum length are not plotted. Would it be reasonable to also let the user specify a minimum length, in addition to the current -Q 'number-of-points' cut off?

I find that if you have the same grid, but different resolution, the -Q 'number of points' may have a big impact.

Pragmatic constructed example:

Plotting and extraction of the -2000 contour from a grid with three different resolutions (see contour_length.pdf)

  • gebco_8m.nc - x/y increment 0.13333
    number of points: N = 121
    Length: 820 km
  • gebco_8m_sampled.nc 1 - x/y increment 0.5
    Number of points: N = 30
    Length: 738 km
  • gebco_8m_sampled.nc 2 - x/y increment 0.05
    Number of points: N = 324
    Length: 822 km

The number of points may be very different, but the length will roughly be the same (at least comparable). Wouldn't length be a nice additional way of specifying contour cut?

contour_length.sh Magnifier (2.8 KB) Andreas, 2017-07-26 23:26

contour_length.pdf (58.2 KB) Andreas, 2017-07-26 23:26

gebco_8m.nc (94.6 KB) Andreas, 2017-07-26 23:26

grdcontour.PNG (38.6 KB) Andreas, 2017-12-23 02:54


#1 Updated by Paul almost 4 years ago

  • Status changed from New to Resolved
  • Assignee set to Paul
  • Target version set to Candidate for next minor release
  • % Done changed from 0 to 100

Finally had time to implement a new -Q in grdcontour and pscontour. Appending a map unit or distance unit switches it to a length test.

#2 Updated by Andreas almost 4 years ago

Tested with a geog grid and worked! However, how would one specify the length when dealing with a projected grid (units in e.g. meter)?

#3 Updated by Paul almost 4 years ago

I revised things so you can give any of the geo units (e, k. f, etc) or c (Cartesian user distances - i.e., whatever the grid units are) and C (projected map distances), whose units depends on PROJ_LENGTH_UNIT. In r19372, see how it goes.

#4 Updated by Andreas almost 4 years ago

I have not been able to test this very thoroughly yet. I tried with a geo grid again;

grdcontour -C1000 -Q100k gebco_8m.nc -D > 100k.gmt

And plotting the resulting lines in QGIS, I do get segments that are < 100 km in length.

#5 Updated by Paul almost 4 years ago

  • Status changed from Resolved to In Progress

Confirmed. Will have a look when I can find time (we are all at AGU).

#6 Updated by Paul over 3 years ago

  • Status changed from In Progress to Resolved

Should be fixed in r19512. I've checked that your case only returns contours >= 100 km.

#7 Updated by Andreas over 3 years ago

Did a quick test. Looks like only contours >= 100 km are returned, but a bunch of valid contours are omitted as well.

The attached screenshot shows a test in QGIS. Left spreadsheet shows all contours and on the right a -Q100k cut-off. The contour with length 148.71 km comes in as the shortest contour above 100 km (on the right), while the spreadsheet on the left (all contours) indicate that there are a lot of other contours longer than 100 km, but less than 148.71 km. (The problem is also valid for contours longer than 148.71 km of course - this is just an example).

Maybe GMT applies both length- and 'number-of-points' cut off?

Hope I'm not spoiling your lutefisk dinner!

#8 Updated by Paul over 3 years ago

If I were having one, I would pay for you to spoil my lutefisk dinner.

Just did the -Q100k and no -Q on your gebco_8m.nc file and examined all contours longer than 100. I got exactly the same 50 using either method. I did

gmt grdcontour -C1000 -Q100k gebco_8m.nc -D > 100k.gmt
gmt grdcontour -C1000 gebco_8m.nc -D > all.gmt
gmt mapproject -G+uk 100k.gmt | gmt convert -El | grep -v '^>' > 100.txt
gmt mapproject -G+uk all.gmt  | gmt convert -El | grep -v '^>' | awk '{if ($4 >= 100) print $0}' > all.txt
diff 100.txt all.txt

and they are identical.

#9 Updated by Andreas over 3 years ago

Oops. I used gmt spatial -Qk to check the distances. I see now that I get a mix of computed lengths and areas when using this. Sorry - thanks Paul.

#10 Updated by Paul over 3 years ago

  • Status changed from Resolved to Closed

Excellent, closing this issue. Merry Christmas, Andreas, thanks for all your hell in 2017, and go easy on the schnapps.

#11 Updated by Andreas over 3 years ago

Thank you - I'll do my best. Merry Christmas, Paul!

Also available in: Atom PDF