Let grdtrend fit model along xx or yy only. #5496

joa-quim · 2021-07-18T15:20:12Z

Sometimes it may be desirable to fit a model that varies only along one direction. An example will be to remove a N-S gradient only. This PR add +x or +y as optional modifiers to -N. However I had some difficulties because the gauss function didn't like to have a set of coords all equal to zero. The solution was to let the last xvar or yvar = 1 and write a function to remove the vertical slice shown in attached image.
Maybe something more clever can be done but meanwhile this now works well.

PaulWessel · 2021-07-19T04:48:10Z

Hi @joa-quim. I agree that finding just x or y trends is very useful. However, I wish you would have first proposed this feature so we could discuss how to do this rather than implementing it, making it harder to suggest changes. I spent a lot of time a few years ago rewriting how trend1d works to allow a wide mix of Fourier and polynomial coefficients, and the plan was to expand this to trend2d and grdtrend. If we add +x and +y then that sets in place things that is hard to undo later (backwards compatibility). Have you tested the new functionality and are there any test scripts so we know it works, etc etc.?

joa-quim · 2021-07-22T16:13:11Z

Yep, trend2d can receive this options as well but didn't look at it. Ofc, the syntax or the technical implementation can be modified but this is how I managed to do it.

Running this script shows it works but I don't know how to turn it into a test that says if 0; 0 then it passed

#!/usr/bin/env bash
#
#	Testing the +x and +y modifiers to grdtrend -N

gmt grdmath -R-1/1/-1/1 -I0.1 X = x.grd
gmt grdtrend x.grd -N3+x -Dlixo_flat_x.grd

gmt grdmath -R-1/1/-1/1 -I0.1 Y = y.grd
gmt grdtrend y.grd -N3+y -Dlixo_flat_y.grd

gmt grdinfo lixo_flat_x.grd -Cn -o4,5 | awk '{print($1 + $2)}'
gmt grdinfo lixo_flat_y.grd -Cn -o4,5 | awk '{print($1 + $2)}'

PaulWessel · 2021-08-23T20:18:58Z

I've decided to accept this PR but will first work on that test script.

PaulWessel · 2021-08-23T20:54:24Z

Wait a bit. I just tried the above script and it fails. it may print 0 0 but the grids are screwed, no? The info says they are both all zeros.

PaulWessel · 2021-08-23T23:08:53Z

Oops, sorry, I see you are using -D to get the residuals. Anyway, I am making some changes so that you dont have to set -N10 just go get 4 terms in y, and improve the docs. I removed the special treatment and just solve for the exact solution with the parameters that are needed only.

PaulWessel

Tested and documented. Good to go, @joa-quim.

PaulWessel · 2021-08-24T00:32:50Z

Since I made many changes I think we may need an independent review on this, @meghanrjones.

maxrjones · 2021-08-24T17:36:23Z

@PaulWessel, is this rounding error within expectations for GMT?

gmt grdmath -R0/1/0/1 -I0.1 X = data.grd
gmt grdtrend data.grd -N2+x -Ddiff.grd -Ttrend.grd

ncdump data.grd

netcdf data {
dimensions:
        x = 11 ;
        y = 11 ;
variables:
        double x(x) ;
                x:long_name = "x" ;
                x:actual_range = 0., 1. ;
        double y(y) ;
                y:long_name = "y" ;
                y:actual_range = 0., 1. ;
        float z(y, x) ;
                z:long_name = "z" ;
                z:_FillValue = NaNf ;
                z:actual_range = 0., 1. ;

// global attributes:
                :Conventions = "CF-1.7" ;
                :title = "Produced by grdmath" ;
                :history = "grdmath -R0/1/0/1 -I0.1 X = data.grd" ;
                :GMT_version = "6.3.0_320bdf7_2021.08.23 [64-bit]" ;
data:

 x = 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 ;

 y = 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 ;

 z =
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 ;
}

ncdump trend.grd

netcdf trend {
dimensions:
        x = 11 ;
        y = 11 ;
variables:
        double x(x) ;
                x:long_name = "x" ;
                x:actual_range = 0., 1. ;
        double y(y) ;
                y:long_name = "y" ;
                y:actual_range = 0., 1. ;
        float z(y, x) ;
                z:long_name = "z" ;
                z:_FillValue = NaNf ;
                z:actual_range = 0., 1. ;

// global attributes:
                :Conventions = "CF-1.7" ;
                :title = "Produced by grdmath" ;
                :history = "grdtrend data.grd -N2+x -Ddiff.grd -Ttrend.grd" ;
                :description = "trend surface" ;
                :GMT_version = "6.3.0_320bdf7_2021.08.23 [64-bit]" ;
data:

 x = 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 ;

 y = 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 ;

 z =
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1,
  0, 0.09999999, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 ;
}

maxrjones · 2021-08-24T17:38:12Z

I would expect these two commands to produce the same trend, since the linear x term should be included in both equations. But, -N3+x results in a trend grid with a constant value:

Example 1:

    gmt grdmath -R0/1/0/1 -I0.1 X = data.grd
    gmt grdtrend data.grd -N2+x -Ddiff.grd -Ttrend.grd

Example 2:

    gmt grdmath -R0/1/0/1 -I0.1 X = data.grd
    gmt grdtrend data.grd -N3+x -Ddiff.grd -Ttrend.grd

maxrjones · 2021-08-24T17:47:42Z

I would expect these two commands to produce the same trend, since the linear x term should be included in both equations. But, -N3+x results in a trend grid with a constant value:

Example 1:
    gmt grdmath -R0/1/0/1 -I0.1 X = data.grd
    gmt grdtrend data.grd -N2+x -Ddiff.grd -Ttrend.grd
Example 2:
    gmt grdmath -R0/1/0/1 -I0.1 X = data.grd
    gmt grdtrend data.grd -N3+x -Ddiff.grd -Ttrend.grd

I misinterpret the documentation for how n_model works when +x|+y is used. Still surprised by the output from the second example.

The cutoff is different for +x or +y since we are doing a linear fit only.

PaulWessel · 2021-08-24T18:44:57Z

I forgot to update the decision on whether we have a trivial solution or not. The check is different now for +x +y. Hopefully this fix is OK.

maxrjones · 2021-08-24T20:01:34Z

Much better. Are the rounding errors acceptable (#5496 (comment))?
There is a problem with -N4 (possibly indexing issue):

gmt grdmath -R0/1/0/1 -I0.1 X X MUL X MUL = data.grd
gmt grdtrend data.grd -N4+x -Ddiff.grd -Ttrend.grd
grdtrend [ERROR]: Option -N: Specify 1-4 model parameters when +x or +y are active

PaulWessel · 2021-08-24T22:04:09Z

Much better. Are the rounding errors acceptable (#5496 (comment))?

Yes. gmt math -Q 0.1 0.09999999 SUB =
1.00000000086e-08

so for 4-byte floats that is as good as it gets.
I fixed the other warning.

maxrjones · 2021-08-25T00:02:40Z

Sorry, but I am skeptical about the least-squares results. For this example with a synthetic grid Z = 0.5 + X + X^2 + X^3, the magnitude of the residuals seems unreasonably large. Am I missing something?

rm data.grd
rm trend.grd
rm diff.grd

gmt begin test png
    gmt grdmath -R0/1/0/1 -I0.1 0.5 = data1.grd
    gmt grdmath -R0/1/0/1 -I0.1 X = data2.grd
    gmt grdmath -R0/1/0/1 -I0.1 X X MUL = data3.grd
    gmt grdmath -R0/1/0/1 -I0.1 X X MUL X MUL = data4.grd
    gmt grdmath data1.grd data2.grd ADD data3.grd ADD data4.grd ADD = data.grd
    gmt grdtrend data.grd -N4+r+x -Ddiff.grd -Ttrend.grd -V
    gmt subplot begin 3x1 -R0/1/0/1 -Ff5c/15c -Bafg -M0.5c
        gmt subplot set 0
            grid="data.grd"
            range=`gmt grdinfo ${grid} -T`
            gmt makecpt -Cwhite,50 ${range} -Z
            gmt grdimage ${grid} -C -JX? -B+t"original"
            gmt colorbar -DJMR+v
        gmt subplot set 1
            grid="trend.grd"
            gmt makecpt -Cwhite,50 ${range} -Z
            gmt grdimage ${grid} -C -JX? -B+t"trend"
            gmt colorbar -DJMR+v
        gmt subplot set 2
            grid="diff.grd"
            range=`gmt grdinfo ${grid} -T`
            gmt makecpt -Cwhite,50 ${range} -Z
            gmt grdimage ${grid} -C -JX? -B+t"diff"
            gmt colorbar -DJMR+v
    gmt subplot end
gmt end show

PaulWessel · 2021-08-25T00:13:17Z

That looks pretty bad. Then I go to master and do the same but fit -N10 which ought to give a good fit and it is just as bad. So something is not working as we expect here. Could you perhaps modify this script to show a crossection at y = 0.5 instead?

maxrjones · 2021-08-25T00:45:30Z

Here you go:

rm data.grd
rm trend.grd
rm diff.grd

gmt begin test png
    gmt grdmath -R0/1/0/1 -I0.1 0.5 = data1.grd
    gmt grdmath -R0/1/0/1 -I0.1 X = data2.grd
    gmt grdmath -R0/1/0/1 -I0.1 X X MUL = data3.grd
    gmt grdmath -R0/1/0/1 -I0.1 X X MUL X MUL = data4.grd
    gmt grdmath data1.grd data2.grd ADD data3.grd ADD data4.grd ADD = data.grd
    gmt grdtrend data.grd -N4+x -Ddiff.grd -Ttrend.grd -V
    gmt grd2xyz -R0/01/0.4/0.5 -o0,2 data.grd > data.txt
    gmt grd2xyz -R0/01/0.4/0.5 -o0,2 trend.grd > trend.txt
    gmt grd2xyz -R0/01/0.4/0.5 -o0,2 diff.grd > diff.txt
    gmt subplot begin 3x1 -Ff5c/15c -Bafg -M0.5c
        gmt subplot set 0
            gmt plot -JX? -B+t"data" -R0/1/0/3.5 -q11: data.txt
        gmt subplot set 1
            gmt plot -JX? -B+t"trend" -R0/1/0/3.5 -q11: trend.txt
        gmt subplot set 2
            gmt plot -JX? -B+t"diff" -R0/1/-1/1 -q11: diff.txt
    gmt subplot end
gmt end show

PaulWessel · 2021-08-25T00:57:18Z

Hm, that is pretty revealing. It is the same that master gives (using -N10 instead of -N4+x).

maxrjones · 2021-08-25T02:22:44Z

Hm, that is pretty revealing. It is the same that master gives (using -N10 instead of -N4+x).

Yes, I think this is not specific to this branch. So this could be merged and the bug fixed as a separate issue.

PaulWessel · 2021-08-25T20:06:58Z

Because this PR (a) adds new capability (that is working) and (b) retains the same problem as documented above as in master, we will merge this PR in and then work on addressing the wider issue in a new PR.

Let grdtrend fit model along xx or yy only.

a94243a

joa-quim requested a review from PaulWessel July 18, 2021 15:20

Merge branch 'master' into grdtrend-x_or_y

b723608

vercel bot temporarily deployed to Preview August 23, 2021 20:20 Inactive

Update docs and add a test script

c40649e

vercel bot temporarily deployed to Preview August 23, 2021 20:55 Inactive

Finalize

3e9cab0

vercel bot temporarily deployed to Preview August 23, 2021 22:46 Inactive

Polish

30e5f51

vercel bot temporarily deployed to Preview August 23, 2021 23:07 Inactive

PaulWessel added the new core module feature PR that implements a new core module feature label Aug 23, 2021

Update grdtrend.rst

6330c67

vercel bot temporarily deployed to Preview August 23, 2021 23:22 Inactive

PaulWessel approved these changes Aug 23, 2021

View reviewed changes

Update grdtrend.rst

320bdf7

vercel bot temporarily deployed to Preview August 24, 2021 03:27 Inactive

maxrjones added the add-changelog Add PR to the changelog label Aug 24, 2021

Update decision on trivial solution

fde18c6

The cutoff is different for +x or +y since we are doing a linear fit only.

vercel bot temporarily deployed to Preview August 24, 2021 18:44 Inactive

Update grdtrend.c

19ef0ea

vercel bot temporarily deployed to Preview August 24, 2021 22:02 Inactive

Merge branch 'master' into grdtrend-x_or_y

184e380

vercel bot temporarily deployed to Preview August 24, 2021 22:09 Inactive

Merge branch 'master' into grdtrend-x_or_y

d809219

vercel bot temporarily deployed to Preview August 25, 2021 20:05 Inactive

PaulWessel merged commit 0b1f96a into master Aug 25, 2021

PaulWessel deleted the grdtrend-x_or_y branch August 25, 2021 20:48

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Let grdtrend fit model along xx or yy only. #5496

Let grdtrend fit model along xx or yy only. #5496

joa-quim commented Jul 18, 2021 •

edited

Loading

PaulWessel commented Jul 19, 2021

joa-quim commented Jul 22, 2021

PaulWessel commented Aug 23, 2021

PaulWessel commented Aug 23, 2021

PaulWessel commented Aug 23, 2021

PaulWessel left a comment

PaulWessel commented Aug 24, 2021

maxrjones commented Aug 24, 2021

maxrjones commented Aug 24, 2021

maxrjones commented Aug 24, 2021

PaulWessel commented Aug 24, 2021

maxrjones commented Aug 24, 2021

PaulWessel commented Aug 24, 2021

maxrjones commented Aug 25, 2021

PaulWessel commented Aug 25, 2021

maxrjones commented Aug 25, 2021

PaulWessel commented Aug 25, 2021

maxrjones commented Aug 25, 2021

PaulWessel commented Aug 25, 2021

Let grdtrend fit model along xx or yy only. #5496

Let grdtrend fit model along xx or yy only. #5496

Conversation

joa-quim commented Jul 18, 2021 • edited Loading

PaulWessel commented Jul 19, 2021

joa-quim commented Jul 22, 2021

PaulWessel commented Aug 23, 2021

PaulWessel commented Aug 23, 2021

PaulWessel commented Aug 23, 2021

PaulWessel left a comment

Choose a reason for hiding this comment

PaulWessel commented Aug 24, 2021

maxrjones commented Aug 24, 2021

maxrjones commented Aug 24, 2021

maxrjones commented Aug 24, 2021

PaulWessel commented Aug 24, 2021

maxrjones commented Aug 24, 2021

PaulWessel commented Aug 24, 2021

maxrjones commented Aug 25, 2021

PaulWessel commented Aug 25, 2021

maxrjones commented Aug 25, 2021

PaulWessel commented Aug 25, 2021

maxrjones commented Aug 25, 2021

PaulWessel commented Aug 25, 2021

joa-quim commented Jul 18, 2021 •

edited

Loading