Hi,
with "r.resamp.stats method=sum", I'm seeing banding. Is this expected
or a case of it including points falling on the northern and eastern
bounds? so the same input cell gets counted for multiple output cells?
i.e. is it doing ( <= x <= , <= y <= ) instead of ( <= x < , <= y < ) ?
#spearfish
g.region rast=elevation.10m
g.region res=131.328
r.resamp.stats in=elevation.10m out=resamp.sum meth=sum
Hamish
Hamish wrote:
with "r.resamp.stats method=sum", I'm seeing banding. Is this expected
or a case of it including points falling on the northern and eastern
bounds? so the same input cell gets counted for multiple output cells?
i.e. is it doing ( <= x <= , <= y <= ) instead of ( <= x < , <= y < ) ?
Unless the destination resolution is an integer multiple of the source
resolution, the total number of source cells used for each destination
cell will vary.
E.g. if the source resolution is 10m, and the destination resolution
is 25m, then the output cells will be based upon grids of 2x2, 2x3,
3x2 and 3x3 source cells, in equal proportions.
For the "sum" aggregate, this will have a significant effect upon the
output, as a cell computed from 3x3 source cells will contain the sum
of over twice as many values as one based upon 2x2 source cells.
The banding effect should go away if you use the -w switch. That will
cause the aggregates to use proportional fractions of each source
cell's value, rather than attributing the entire source cell to
whichever destination cell contains the source cell's centre. E.g. in
the 10m->25m case, each destination cell will use values from 2.5x2.5
source cells (i.e. it will use 2x2 values with a weight of 1, 2x1+1x2
values with a weight of 0.5, and 1x1 value with a weight of 0.25).
Using the -w switch will probably also be a lot slower.
If the banding occurs with integer scale factors, or with -w, let me
know.
--
Glynn Clements <glynn@gclements.plus.com>
Glynn Clements wrote:
Hamish wrote:
> with "r.resamp.stats method=sum", I'm seeing banding. Is this
> expected or a case of it including points falling on the northern
> and eastern bounds? so the same input cell gets counted for multiple
> output cells? i.e. is it doing ( <= x <= , <= y <= ) instead of ( <=
> x < , <= y < ) ?Unless the destination resolution is an integer multiple of the source
resolution, the total number of source cells used for each destination
cell will vary.E.g. if the source resolution is 10m, and the destination resolution
is 25m, then the output cells will be based upon grids of 2x2, 2x3,
3x2 and 3x3 source cells, in equal proportions.
..
The banding effect should go away if you use the -w switch. That will
cause the aggregates to use proportional fractions of each source
cell's value, rather than attributing the entire source cell to
whichever destination cell contains the source cell's centre.
Yes, adding the -w flag did the trick now it looks good. Both slope and
curvature derivatives are perfectly clean as well.
thanks,
Hamish
ps - I've added metadata (hist/) support to r.resamp.stats & .interp,
and they now copy over the source map's color table, for all cases but
method=sum.
On Monday 30 October 2006 20:00, Hamish wrote:
Glynn Clements wrote:
> Hamish wrote:
> > with "r.resamp.stats method=sum", I'm seeing banding. Is this
> > expected or a case of it including points falling on the northern
> > and eastern bounds? so the same input cell gets counted for multiple
> > output cells? i.e. is it doing ( <= x <= , <= y <= ) instead of ( <=
> > x < , <= y < ) ?
>
> Unless the destination resolution is an integer multiple of the source
> resolution, the total number of source cells used for each destination
> cell will vary.
>
> E.g. if the source resolution is 10m, and the destination resolution
> is 25m, then the output cells will be based upon grids of 2x2, 2x3,
> 3x2 and 3x3 source cells, in equal proportions...
> The banding effect should go away if you use the -w switch. That will
> cause the aggregates to use proportional fractions of each source
> cell's value, rather than attributing the entire source cell to
> whichever destination cell contains the source cell's centre.Yes, adding the -w flag did the trick now it looks good. Both slope and
curvature derivatives are perfectly clean as well.thanks,
Hamishps - I've added metadata (hist/) support to r.resamp.stats & .interp,
and they now copy over the source map's color table, for all cases but
method=sum.
Great job Gylnn on implementing this module. I lost track of its progress
while writting up a recent proposal for a very similar topic. I can confirm
the behavior that Hamish has reported. I'll post back with a comparison
between the output from the new module, and the vector/raster hack that I had
used in the past: http://casoilresource.lawr.ucdavis.edu/drupal/node/275
Cheers,
--
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341
Hi,
I've added r.resamp.stats and r.resamp.interp to the GUI menus.
(Raster->Devel)
should r.resample be tagged for future removal? (same result as with:
'r.resamp.interp method=nearest' ???).
Hamish
Hamish wrote:
I've added r.resamp.stats and r.resamp.interp to the GUI menus.
(Raster->Devel)should r.resample be tagged for future removal? (same result as with:
'r.resamp.interp method=nearest' ???).
No.
r.resample uses libgis' built-in resampling, so it is guaranteed to
resample in an identical manner to other commands.
"r.resamp.interp method=nearest" performs the resampling itself, so it
could produce slightly different results due to rounding errors. The
method=nearest option was only included for completeness; if you want
nearest-neighbour resampling, r.resample is preferred.
--
Glynn Clements <glynn@gclements.plus.com>
Glynn Clements wrote on 11/02/2006 05:26 AM:
Hamish wrote:
I've added r.resamp.stats and r.resamp.interp to the GUI menus.
(Raster->Devel)should r.resample be tagged for future removal? (same result as with:
'r.resamp.interp method=nearest' ???).
No.r.resample uses libgis' built-in resampling, so it is guaranteed to
resample in an identical manner to other commands."r.resamp.interp method=nearest" performs the resampling itself, so it
could produce slightly different results due to rounding errors. The
method=nearest option was only included for completeness; if you want
nearest-neighbour resampling, r.resample is preferred.
Honestly, this sounds confusing to me. Having NN method twice which produces
slightly different results, is hard to explain to users.
Markus
Markus Neteler wrote:
>> I've added r.resamp.stats and r.resamp.interp to the GUI menus.
>> (Raster->Devel)
>>
>> should r.resample be tagged for future removal? (same result as with:
>> 'r.resamp.interp method=nearest' ???).
>>
>
> No.
>
> r.resample uses libgis' built-in resampling, so it is guaranteed to
> resample in an identical manner to other commands.
>
> "r.resamp.interp method=nearest" performs the resampling itself, so it
> could produce slightly different results due to rounding errors. The
> method=nearest option was only included for completeness; if you want
> nearest-neighbour resampling, r.resample is preferred.Honestly, this sounds confusing to me. Having NN method twice which produces
slightly different results, is hard to explain to users.
It's largely inevitable.
Even if r.resample is removed, it can (mostly) be duplicated with e.g.
"r.mapcalc dst = src", which also isn't guaranteed to produce the same
results as r.resamp.interp. Similarly for any other module which can
be made to perform a "copy" operation by some setting of its
parameters (r.mfilter with a 1x1 kernel, maybe?).
I could remove method=nearest from r.resamp.interp, but if it's ever
extended to support more advanced interpolation schemes with
parameters, the same issue could arise.
Ultimately, nearest-neighbour resampling will always have boundary
cases where the centre of an output cell falls on the boundary between
two or four input cells. In that situation, the end result will depend
upon everything which can affect the rounding error: the algorithm
used, the compiler, the compilation switches, the architecture, etc.
--
Glynn Clements <glynn@gclements.plus.com>
On Thursday 02 November 2006 17:03, Glynn Clements wrote:
Markus Neteler wrote:
> >> I've added r.resamp.stats and r.resamp.interp to the GUI menus.
> >> (Raster->Devel)
> >>
> >> should r.resample be tagged for future removal? (same result as with:
> >> 'r.resamp.interp method=nearest' ???).
> >
> > No.
> >
> > r.resample uses libgis' built-in resampling, so it is guaranteed to
> > resample in an identical manner to other commands.
> >
> > "r.resamp.interp method=nearest" performs the resampling itself, so it
> > could produce slightly different results due to rounding errors. The
> > method=nearest option was only included for completeness; if you want
> > nearest-neighbour resampling, r.resample is preferred.
>
> Honestly, this sounds confusing to me. Having NN method twice which
> produces slightly different results, is hard to explain to users.It's largely inevitable.
Even if r.resample is removed, it can (mostly) be duplicated with e.g.
"r.mapcalc dst = src", which also isn't guaranteed to produce the same
results as r.resamp.interp. Similarly for any other module which can
be made to perform a "copy" operation by some setting of its
parameters (r.mfilter with a 1x1 kernel, maybe?).I could remove method=nearest from r.resamp.interp, but if it's ever
extended to support more advanced interpolation schemes with
parameters, the same issue could arise.Ultimately, nearest-neighbour resampling will always have boundary
cases where the centre of an output cell falls on the boundary between
two or four input cells. In that situation, the end result will depend
upon everything which can affect the rounding error: the algorithm
used, the compiler, the compilation switches, the architecture, etc.
This is deviating someone from the NN issues... but:
I have been noticing some odd edge effects when using r.resamp.stats
method=mode .... i.e. i am getting a "ring" of cells with 0-value around the
actual data. Within the "ring" the results appear to match those obtained
with an alternate raster-vector-starspan-raster approach to doing the same
operation.
a graphical summary can be found here:
http://169.237.35.250/~dylan/temp/method_mode-test.png
Details on the raster-vector-raster approach are listed here:
http://casoilresource.lawr.ucdavis.edu/drupal/node/275
cheers,
--
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341
Dylan Beaudette wrote:
I have been noticing some odd edge effects when using r.resamp.stats
method=mode .... i.e. i am getting a "ring" of cells with 0-value around the
actual data. Within the "ring" the results appear to match those obtained
with an alternate raster-vector-starspan-raster approach to doing the same
operation.a graphical summary can be found here:
http://169.237.35.250/~dylan/temp/method_mode-test.png
I wasn't making the source window large enough, with the result that
edge cells were using garbage values from outside of the source
buffer.
[r.resamp.stats was derived from r.resamp.interp. The window-setting
code was left untouched, but r.resamp.interp only needs source values
corresponding to the destination cell's centre +/- 2 source cells,
while r.resamp.stats needs source values for the whole of the
destination cell's area.]
I've committed the attached patch to CVS. Let me know if you still
encounter problems.
--
Glynn Clements <glynn@gclements.plus.com>
r.resamp.stats.diff (1.15 KB)
On Friday 03 November 2006 14:51, Glynn Clements wrote:
Dylan Beaudette wrote:
> I have been noticing some odd edge effects when using r.resamp.stats
> method=mode .... i.e. i am getting a "ring" of cells with 0-value around
> the actual data. Within the "ring" the results appear to match those
> obtained with an alternate raster-vector-starspan-raster approach to
> doing the same operation.
>
> a graphical summary can be found here:
> http://169.237.35.250/~dylan/temp/method_mode-test.pngI wasn't making the source window large enough, with the result that
edge cells were using garbage values from outside of the source
buffer.[r.resamp.stats was derived from r.resamp.interp. The window-setting
code was left untouched, but r.resamp.interp only needs source values
corresponding to the destination cell's centre +/- 2 source cells,
while r.resamp.stats needs source values for the whole of the
destination cell's area.]I've committed the attached patch to CVS. Let me know if you still
encounter problems.
Glynn, others:
The latest patch looks like a winner!
I have tested out the latest r.resamp.stats, comparing the results betwee:
1. starspan aggregation and r.resamp.stats
2. r.resamp.stats and default NN resampling
I have posted a detailed set of diagnostics here:
http://casoilresource.lawr.ucdavis.edu/drupal/node/275
It looks like r.resamp.stats is returning good results, which are very close
to those returned by the starspan hack I had used before.
Cheers,
--
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341
I asked earlier if r.resample should be retired, the answer was
"no, keep it" as `r.resamp.interp method=nearest` isn't guaranteed to
give exactly the same result. (but `r.mapcalc new=old` will match r.resample?)
how about this one:
Should r.bilinear be tagged for retirement after GRASS 6.x?
(due to `r.resamp.interp method=bilinear`)
Hamish
Hamish wrote:
I asked earlier if r.resample should be retired, the answer was
"no, keep it" as `r.resamp.interp method=nearest` isn't guaranteed to
give exactly the same result. (but `r.mapcalc new=old` will match r.resample?)
Just about any r.* command which has a "copy mode" (i.e. some set of
arguments which results in the output being equal to the input) will
behave like r.resample. The core of r.resample is:
for (row = 0; row < nrows; row++)
{
if (verbose)
G_percent (row, nrows, 2);
if (G_get_raster_row (infd, rast, row, data_type) < 0)
exit(1);
if (G_put_raster_row (outfd, rast, out_type) < 0)
exit(1);
G_mark_raster_cats (rast, ncols, &cats, data_type);
}
[Literal cut-and-paste from raster/r.resample/main.c]
IOW, it's a get-row/put-row loop, which relies upon the automatic
resampling performed by G_get_raster_row() etc. Any GRASS module which
reads a raster map using the current region will have the map
resampled in exactly the same way.
how about this one:
Should r.bilinear be tagged for retirement after GRASS 6.x?
(due to `r.resamp.interp method=bilinear`)
Yes.
Or replaced by a script which prints a "this module is deprecated"
warning then runs "r.resamp.interp method=bilinear ...".
--
Glynn Clements <glynn@gclements.plus.com>
On Wednesday 08 November 2006 23:38, Glynn Clements wrote:
Hamish wrote:
> I asked earlier if r.resample should be retired, the answer was
> "no, keep it" as `r.resamp.interp method=nearest` isn't guaranteed to
> give exactly the same result. (but `r.mapcalc new=old` will match
> r.resample?)Just about any r.* command which has a "copy mode" (i.e. some set of
arguments which results in the output being equal to the input) will
behave like r.resample. The core of r.resample is:for (row = 0; row < nrows; row++)
{
if (verbose)
G_percent (row, nrows, 2);
if (G_get_raster_row (infd, rast, row, data_type) < 0)
exit(1);
if (G_put_raster_row (outfd, rast, out_type) < 0)
exit(1);
G_mark_raster_cats (rast, ncols, &cats, data_type);
}[Literal cut-and-paste from raster/r.resample/main.c]
IOW, it's a get-row/put-row loop, which relies upon the automatic
resampling performed by G_get_raster_row() etc. Any GRASS module which
reads a raster map using the current region will have the map
resampled in exactly the same way.> how about this one:
>
> Should r.bilinear be tagged for retirement after GRASS 6.x?
> (due to `r.resamp.interp method=bilinear`)Yes.
Or replaced by a script which prints a "this module is deprecated"
warning then runs "r.resamp.interp method=bilinear ...".
I think that this last option is a good migration path to a more coherent set
of resampling modules.
Thanks again Glynn for setting these up. However, I have noticed odd behavior
with r.resample.stats when the native resolution == the current resolution:
i.e. there is no difference in start/final cell size. I will start a new
thread on this topic.
Cheers,
--
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341
Glynn Clements wrote:
> Should r.bilinear be tagged for retirement after GRASS 6.x?
> (due to `r.resamp.interp method=bilinear`)Yes.
Or replaced by a script which prints a "this module is deprecated"
warning then runs "r.resamp.interp method=bilinear ...".
warning added.
note r.bilinear has these two options which r.resamp.interp doesn't:
north specific input value to be assigned to the north and/or south poles for longitude-latitude grids
east specific input value to be assigned to the north and/or south poles for longitude-latitude grids
Hamish
Hamish wrote:
> > Should r.bilinear be tagged for retirement after GRASS 6.x?
> > (due to `r.resamp.interp method=bilinear`)
>
> Yes.
>
> Or replaced by a script which prints a "this module is deprecated"
> warning then runs "r.resamp.interp method=bilinear ...".warning added.
note r.bilinear has these two options which r.resamp.interp doesn't:
north specific input value to be assigned to the north and/or south poles for longitude-latitude grids
east specific input value to be assigned to the north and/or south poles for longitude-latitude grids
Hmm. I suspect that the region enlargement will cause r.resamp.interp
will fail with "Illegal latitude for North" (or South) if you try to
resample a lat/lon map which includes either pole.
I'm starting to wonder whether we would be better off without the
latitude check, leaving it to up the user instead.
--
Glynn Clements <glynn@gclements.plus.com>
On Monday 13 November 2006 00:36, Glynn Clements wrote:
Hamish wrote:
> > > Should r.bilinear be tagged for retirement after GRASS 6.x?
> > > (due to `r.resamp.interp method=bilinear`)
> >
> > Yes.
> >
> > Or replaced by a script which prints a "this module is deprecated"
> > warning then runs "r.resamp.interp method=bilinear ...".
>
> warning added.
>
>
> note r.bilinear has these two options which r.resamp.interp doesn't:
>
> north specific input value to be assigned to the north and/or south
> poles for longitude-latitude grids east specific input value to be
> assigned to the north and/or south poles for longitude-latitude gridsHmm. I suspect that the region enlargement will cause r.resamp.interp
will fail with "Illegal latitude for North" (or South) if you try to
resample a lat/lon map which includes either pole.I'm starting to wonder whether we would be better off without the
latitude check, leaving it to up the user instead.
Would a note of caution in the man pages do the trick ?
--
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341
Glynn Clements wrote:
Hmm. I suspect that the region enlargement will cause r.resamp.interp
will fail with "Illegal latitude for North" (or South) if you try to
resample a lat/lon map which includes either pole.I'm starting to wonder whether we would be better off without the
latitude check, leaving it to up the user instead.
I think the latitude check is quite valuable. I don't see a reason to
allow incorrect results just because the user asked to do something
impossible (probably by accident).
In the case of r.resamp.interp it would be nice if the module could
silently convert affected cells to NULL instead of spitting out error
messages of course.
Hamish
