Bug #618

grdfilter boxcar with -distance returns zero

Added by Joseph almost 3 years ago. Updated almost 3 years ago.

Status:ClosedStart date:2014-09-09
Priority:NormalDue date:
Assignee:Paul% Done:

100%

Category:-
Target version:Candidate for next bugfix release
Affected version:5.1.1 Platform:

Description

grdfilter -Fb-600 returns a grid of zero.

grdfilt -Fm-600 works, as does -Fg-600 and so on.

I'll speculate the bug is that -B-600 is subtracting the mean from the mean rather than the mean from the value.

Since median is ofter about the same as mean value, that is the work around

History

#1 Updated by Paul almost 3 years ago

  • Status changed from New to In Progress

I cannot reproduce your zero grid with -Fb-600. I tried this:
grdmath -R0/100/-30/60 -I1 -fg 3 3 NRAND = t.nc [some random numbers centered on 30 with sigma = 3]
grdfilter -Fb-600 t.nc -Gr.nc -V -D3 [filter with mean and return residuals]

Find that residuals are now centered around 0 as the mean was removed. Do you get zero on this or is there another problem?

grdinfo r.nc
r.nc: Title: t.nc
r.nc: Command: grdfilter -V -Fb-600 t.nc -Gr.nc -D3
r.nc: Remark:
r.nc: Gridline node registration used
r.nc: Grid file format: nf (# 18) GMT netCDF format (float) (COARDS-compliant) [DEFAULT]
r.nc: x_min: 0 x_max: 100 x_inc: 1 name: longitude [degrees_east] nx: 101
r.nc: y_min: -30 y_max: 60 y_inc: 1 name: latitude [degrees_north] ny: 91
r.nc: z_min: -12.2926673889 z_max: 11.8014545441 name: z
r.nc: scale_factor: 1 add_offset: 0

#2 Updated by Joseph almost 3 years ago

hmmm, your example works, It even works on my data. OK I found the problem: small windows..

If filter width is less than 1 "pixel", aka 111km on a 1 deg grid, the output is as expected 0.

With low pass filters, ( aka positive distance), it works with even a tiny 3x3 convolution, (filter widths between 1 and 2 pixels).

However with a negative distance it needs at least 5x5.

I am not concerned that quintiles of 9 points are weird, and I can accept 5x5 min, if that is intended behavior.

Perhaps add to documentation...

#3 Updated by Paul almost 3 years ago

The reported 3x3 does not mean it is a 3x3 filter; it means the max arrays for filter weights are 3x3. In the cases you get 0 it means you only averaged a single point and hence the mean = that point and residuals would be zero.

#4 Updated by Joseph almost 3 years ago

I agree with "averaged a single point and hence the mean = that point and residuals would be zero. , but the " MAX" in "it means the >>MAX<< arrays for filter weights are 3x3" is a big deal.

I was thinking that 3x3, 7x7, or whatever filter size >>displayed<< by the grdfilter output, was what it >>used<<, and hence it was grabbing a matching window out of the data.

If it doesn't do that, don't all the convolution filters have interesting side effects on geographic grids?

Nothing exactly wrong with that, but again even a brief mention of that on the man page would be very useful; e,g, "If the filter wight matrix is bigger than the window, it extends the window with NaN", or whatever it does.

In the big picture, I suppose people are expecting something loosely called "auto-magic equal area" filtering, with the usual problems of projecting, and this perhaps involved. Speaking of that, a 600km filter window at high latitudes is a whole lot of pixels, does the filter matrix get bigger, or what happens there?

#5 Updated by Paul almost 3 years ago

  • Status changed from In Progress to Resolved
  • Assignee set to Paul
  • Target version set to Candidate for next bugfix release
  • % Done changed from 0 to 90

Of course grdfilter is dealing with the change in the number of nodes that fit inside a 600 km circle on the Earth as a function of latitude. If your grid covers the pole you will see that max nx filter width equals your grid width. The best we can do is to modify this message - it is not really that useful to users and was there for us developers to see that indeed the filter nx ny rectangle got set correctly. Because of making sure we are not too small we use ceil and floor and ring for these things so while it may say 3x3 it may indeed be the case that only the center point is within the actual radius you gave (certainly the 4 corners of the 3x3 would ALWAYS be outside the circle, so it is a question of whether the 4 other points are inside r or not. I will improve the verbose messages.

#6 Updated by Paul almost 3 years ago

  • Status changed from Resolved to Closed
  • % Done changed from 90 to 100

I have moved the specific filter nx,ny message for developers to long verbose. In r13517. Closing this issue.

Also available in: Atom PDF