Add grdclip functionality to grdmath

[anbj]

Why is the functionality in grdclip not included (or exclusively) in grdmath?

Today I had to do some calculations with grdmath and then set all cell values < 0 to NaN. To do this I had to:

grdmath ........
grdclip -Sb0/NAN .........

Ideally, this should be possible to do with a single call to grdmath.
In my example only -Sb from grdclip was used, but all the options in grdclip could be included in grdmath.

[PaulWessel]

Because grdclip predates grdmath. But I agree we could add those features to grdmath via new operators, e.g.,

low replace CLIPL (-Sb)
high replace CLIPH (-Sa)
low high replace RANGEFILL (-Si)
value replace REPLACE (-Sr)

Since for backwards reasons we cannot remove grdclip we now would have to keep grdclip around but perhaps

1. Not advertise it (remove from docs and examples)
2. Replace it with short program that calls grdmath instead.

[anbj]

That explains it; I was wondering why grdclip existed - all of this could(/should) be done in grdmath.

I think your suggestions are good.

[joa-quim]

grdmath is an impossible module to wrap. But no problems making grdclip a mock-up module. One of these days I'll have to implement it directly in Julia, which will be much more efficient (no need to care with the pad and no memory copies).

[PaulWessel]

I think pyGMT also are not able to wrap grdmath and gmtmath. But we did so in GMT/MEX but of course there we have the text command to work with.

[joa-quim]

Wrapping is not a problem (Julia has it). The problem is to give it a friendly syntax.

[Esteban82]

Just a question. With grdmath it would be possible to reclassify several ranges of values like this:

gmt grdclip classes.nc -Gnew_classes.nc -Si25/30/99 -Si35/39/55 -Sr17/11 -Sb10/0 -V

[PaulWessel]

That is correct, but that is what grdclip allows as well (-Sr and -Si are repeatable options in grdclip).

[PaulWessel]

Also, @anbj, note you could do what you wanted in grdmath in a different way:

grdmath ..... your calculations ... DUP 0 GE 1 NAN MUL = result.grd

Duplicate the stack, compare to 0 and get another grid that is 0 (false) for z < 0 and 1 (true) for z >= 0, then turn those 0s into NaNs, and multiply the 1 and NaN grid with your data, which turns things into what they were or NaN.

[Esteban82]

Thanks.

BTW, I had made a proposal (#4573) to improve grdclip.

[PaulWessel]

Thanks for the reminder. It remains a tricky one and isn't super general. if your break values were computable it might be simpler (meaning clean multiples in the ranges rather than 3, 12, 22 etc. Let's say you want the slope values to be in three groups: 0-10,10-20,20-30 and azimuths in 8 groups.

gmt grdmath slopes.grd 10 DIV FLOOR 10 MUL = slope_in_10s.grd
gmt grdmath azimuth.grd 45 DIV RINT 45 MUL = azim_in_45.grd

But you are still up against the basic idea that this is a 2-D CPT where both azim and slope needs to be looked up from a 2-D CPT. Since we have no such thing it because a lot of scripting, like @joa-quim said.

I imagine the simplest solution for the end user would a module like this

grdvectorimage gridx gridy [-A -R -J -B ..] -D2DCPT.cpt2 -pdf map

where -A is like in grdvector (so take slope and azim) and -D is the beast of a 2-D CPT. In creating the image we need both r and a to determine the rgb from your 2-D cpt somehow, and we may also allow -Iintens.grd. So most of this is like grdimage except

1. 2 input data grids
2. A 2-D CPT

Looking at your color wheel, seems the hue is set via azimuth and then it is saturation changes with slope. Could you simulate this with turning your slope into the appropriate intensities and use grdimage?

[anbj]

Also, @anbj, note you could do what you wanted in grdmath in a different way:

grdmath ..... your calculations ... DUP 0 GE 1 NAN MUL = result.grd

Duplicate the stack, compare to 0 and get another grid that is 0 (false) for z < 0 and 1 (true) for z >= 0, then turn those 0s into NaNs, and multiply the 1 and NaN grid with your data, which turns things into what they were or NaN.

Thanks Paul, for this enlightening tour de force of RPN magic. Unfortunately, I never had an HP calculator.

[joa-quim]

1. Not advertise it (remove from docs and examples)

No, this would be very confusing to the people that know it exists and would not find the manual anymore. And besides, grdclip manual is much simpler to consult than the one for grdmath.

[Esteban82]

I agree with Joaquim. I also think that the functions of grdclip within grdmath will be lost and we may have users asking for them.

[anbj]

I understand that backward compatibility is/may be/is chosen to be important; but doesn't it create a plethora of technical debt? (That applies to all other work that is done in gmt as well - I'm surprised to see how often an 'old way of doing things' is preserved as an option, i.e. backwards compatibility). If people miss it, they need to consult the changelog and find out what happened.

[maxrjones]

@weiji14 created a proof of concept for a PyGMT wrapper of grdmath in GenericMappingTools/pygmt#1527.

I agree with @joa-quim about continuing to advertise grdclip, even if internally it is changed to a shortcut for new grdmath functionality. RPN is not easy for lots of people whereas grdclip is relatively simple.

[anbj]

@PaulWessel, what do you think? Closing time? Let's get rid of it if it's not worth doing.

[PaulWessel]

Yes, please do.