Dear All,
I have a global, geographic grid (regularly spaced, pixel-registered, 0.5-degree data) in which I would like to calculate a spatial moving statistic, where each pixel in an output grid (with the same global extent and 0.5 degree spacing as the input grid) is equal to the sum of all of the pixels within a 1500 km radius of the center of each pixel in the input grid.
This seems like a job for grdfilter but I am not quite sure how to implement it. Reading the documentation for grdfilter, I’m not sure what would be the appropriate -Ddistance_flag.
Could anyone give me an example? I would be so grateful!
Many thanks for the help,
Jed
The flag is easy. If grid is geog and you want all nodes inside a certain radius than -D4 is the more accurate (and slow) but -D3 could be good enough (and faster).
However, grdfilter doesn’t have a option to calculate the sum of all cells. The closes is the Boxcar (mean). But I suspect you would be surprised with the result of the sum of all cells inside a neighborhood of 1500 km. You see, further north you go, more cells fit inside a circle. Earth being round keeps playing tricks on us.
grdfilter -D4 will correctly assign variable weight depending on latitude. So summing nodes without that weight makes no sense since nodes at 60N north would count twice as much as nodes on the equator. I need to explore this a bit more but we might be able to add a SUM operator to grdfilter.
Thanks for the replies Joaquim and Paul!
Indeed, I can get close to what I want, but grdfilter does not have a SUM operator, and neither does any other GMT tool that does what I would like it to do, that I could think of.
But if a boxcar filter calculates a mean, then it must also calculate a sum, no? I would be so grateful if there might be an easy fix to add a SUM operator to grdfilter. I looked at the code myself but am not much of a C programmer and wasn’t really sure where to begin.
Thanks again for your replies and help!
Jed
I don’t understand what kind of operator will that be. A kind of weighted sum? For what?
Well, it would be a weighted sum of the nodes inside, where the weights are related to the areas. We do this already, then multiply those weights with the filter weights. That is the weight we multiply the node value with and sum up (sum_wz) as well as sum up separately (sum_w). The output is sum_wz/sum_w but if the filter was a sum filter then it would be a boxcar filter without that final division I think.
Yes, what I want is what Paul suggests: just sum_wz without the final division.
I don’t need a weighted sum, I just want to get the sum of all of the values – in my case it is gridded dataset of total lightning strokes – within a 1500km radius of each grid node. The quantity I want to sum should not be multiplied by a weight, I want the total to just be the gross number of strokes (not per unit area).
Hm, that means for the same number of lighting strokes per unit area you will end up with much higher counts as a function of latitude. I don’t think that is a good idea. GMT won’t implement something that is wrong. You could always force the grid to be seen as Cartesian with -f0-1f and then we operate on Cartesian and all nodes have equal weight.
I think I have figured this out now. Thank you both for your help and advice. Cheers, Jed