Gradient caluclation with grdgradient

Added by Lukas 6 months ago

Hi,

I run into a minor issue calculating the gradient of a certain field (e.g. velocity on a grid). To test it, I tried to use a fairly simple function (cosine in X) and calculate the gradient in x -direction (i.e. -A90). The pattern looks ok (i.e. a sinus function), however, the amplitude is not correct (should go from -1 to 1). I attached the file and added the script below.

Further on, this seems to be an even bigger issue, when trying to apply it to a velocity field on a sphere (-R0/359.5/89.75/89.75). In this case the gradient fails entirely.

Any suggestions regarding this issue? I would be very thankful for any comments.

Thank you very much.

Best,

Lukas

—--—--—-
#!/bin/bash

#reg=-R0/359.5/-89.75/89.75
#inc=-I.5/.5
#period=180.0

reg=-R0/8.0/0/8.0
inc=-I.05/.05
period=4.0

direction=${1-90}

rm *.ps

grdmath $reg $inc 2 PI MUL $period DIV X MUL COS = output.grd

  1. Get colorscale
    grd2cpt -E21 -Chaxby output.grd > cmp.cpt
  2. projection
    proj=-JX3
    #proj=-JN180/5i
    #proj=-JQ180/5i
    #ann=-Ba100f50:x": #WesN
    ann=-Ba60f10:x":

grdgradient output.grd -A$direction -Ggradient.grd

grdimage output.grd $proj $reg -Ccmp.cpt -K > output.ps
psscale -Ccmp.cpt -O -K -D1.5/-1/3.5/.2h -B.2 >> output.ps
psbasemap -O $reg $proj $ann >> output.ps
modifybb output.ps

grd2cpt -E21 -Chaxby gradient.grd >cmp2.cpt
grdimage gradient.grd $proj $reg -Ccmp2.cpt -K > gradient.ps
psscale -Ccmp2.cpt -O -K -D1.5/-1/3.5/.2h -B.2 >> gradient.ps
psbasemap -O $reg $proj $ann >> gradient.ps
modifybb gradient.ps

#grd2map gradient.grd
#grd2map output.grd

\epsmerge -par -x 2 -y 1 --orientation Landscape --print --postscript \
output.ps gradient.ps > final.ps 2> /dev/null

gv final.ps

test_gradient (1 KB)


Replies (5)

RE: Gradient caluclation with grdgradient - Added by Paul 6 months ago

I disagree. The derivative of cos(2*pi*x/L) is -(2*pi/L)*sin(2*pi*x/L) so there is a factor of 2*pi/L ~ 1.57 which you are seeing.

RE: Gradient caluclation with grdgradient - Added by Lukas 6 months ago

Oh, yes. You are right. That only leaves the issue when applying it to a sphere.

RE: Gradient caluclation with grdgradient - Added by Paul 6 months ago

Give it a shot. We do consider the latitude when doing these things, but of course these are finite differences modified by shrinking dx as we go towards the poles. Important that your grid is COARDS compliant and identifies itself as geographic (i.e., grdinfo will report degrees_east etc.)

RE: Gradient caluclation with grdgradient - Added by Lukas 6 months ago

Ok, that seems to make sense.

As I haven't worked on a sphere before, I am just curious if the result is correct (because it doesn't look right).

I added the script for a field on as sphere, which varies with a cosine function only in longitude by:

grdmath -R0/359.5/89.75/89.75 -I.5/.5 2 PI MUL 180 DIV X MUL COS = output.grd
Which should identify it as a geographic field.

So taking the gradient along the 90 degrees azimuth I expected it to look similar.

It works at least for the pattern using:
grdgradient output.grd -A90 -Lg -Nt -Ggradient.grd

However, considering the -M option it looks entirely different, with variations at the pol and no gradient along the equator:
grdgradient output.grd -A$direction -Lg -M -Ggradient.grd

Is this what I should expect?


#!/bin/bash

reg=-R0/359.5/-89.75/89.75
inc=-I.5/.5
period=180.0

direction=${1-90}

rm *.ps

grdmath $reg $inc 2 PI MUL $period DIV X MUL COS = output.grd

  1. Get colorscale
    grd2cpt -E21 -Chaxby output.grd > cmp.cpt
    proj=-JQ180/5i
    ann=-Ba60f10:x":

grdgradient output.grd -A$direction -Lg -Nt -Ggradient.grd

grdimage output.grd $proj $reg -fg -Ccmp.cpt -K > output.ps
psscale -Ccmp.cpt -O -K -D1.5/-1/3.5/.2h -B.2 >> output.ps
psbasemap -O $reg $proj $ann >> output.ps

grd2cpt -E21 -Chaxby gradient.grd >cmp2.cpt
grdimage gradient.grd $proj $reg -fg -Ccmp2.cpt -K > gradient.ps
psscale -Ccmp2.cpt -O -K -D1.5/-1/3.5/.2h -B.1 >> gradient.ps
psbasemap -O $reg $proj $ann >> gradient.ps

RE: Gradient caluclation with grdgradient - Added by Paul 6 months ago

Remember GMT is not just for geographic data so giving just -R0/359.5/-89.75/89.75 is not enough information to know these are degrees. I think you mean -Rg -I30m -r, for instance, but you can stick with what you used provided you add -fg (force geographic) and I am guessing you want to -r for pixel registration.

(1-5/5)