[GRASS-user] r.surf.nnbathy interpolation to whole region

GRASS GRASS user
1

Dear list,

I have successfully included natural neighbor interpolation in my mapping routine via r.surf.nnbathy: it seems slower (is that OK?) than v.surf.idw but sounds a lot more robust from what I have read around -- this is very useful if the interpolated raster is to be used for further analysis, beyond visualization.

Currently, the interpolated raster output I get from r.surf.nnbathy is confined to the convex hull delimited by the outermost input points/cells. Would there be a way to extend the interpolated raster to the whole current region as in v.surf.idw (in principle this seems possible: outer Voronoi cells just get bigger)?

By looking at r.surf.nnbathy code and nnbathy --help, it seems that nnbathy should produce a rectangular grid of <cols> by <rows>. What am I missing?

Thanks and regards,

Luigi

2

Luigi Ponti wrote:

I have successfully included natural neighbor interpolation in my
mapping routine via r.surf.nnbathy: it seems slower (is that OK?)
than v.surf.idw but sounds a lot more robust from what I have read
around

I'd only like to stress that it really depends on what your input data
looks like and what the purpose is. Each interpolation method has it's
cons and pros. There is no single most robust method for all applications.

-- this is very useful if the interpolated raster is to be
used for further analysis, beyond visualization.

Currently, the interpolated raster output I get from r.surf.nnbathy
is confined to the convex hull delimited by the outermost input
points/cells. Would there be a way to extend the interpolated raster
to the whole current region as in v.surf.idw

1. Note the hardcoded "-W 0" in the script. You can set it to something
smaller if you want to extrapolate beyond the convex hull. However, I
don't think this makes much sense - set "-W -2" and see how bad the
output DEM looks on the edges. And it can't look better, for how the
triangulation works. Triangulation is restricted to the convex hull
made of input points.

2. Even if your input raster is bigger than your current region,
r.stats, which is used in r.surf.nnbathy script to create the x,y,z
input for nnbathy, discards input raster's cells outside the current
region. Thus, the convex hull is created only for non-null input points
falling into the current region. And the convex hull shape will depend
on how the input non-null cells falling into the current region are
distributed. If you want a perfectly rectangular convex hull, make sure
there is a non-null input raster cell in each region's corner.

I'll think of automating this. Next vacation maybe :).

(in principle this
seems possible: outer Voronoi cells just get bigger)?

I don't exactly understand what you mean. Could you elaborate please?

By looking at r.surf.nnbathy code and nnbathy --help, it seems that
nnbathy should produce a rectangular grid of <cols> by <rows>. What
am I missing?

Yes, it does produce a rectangular grid. Only that all the cells
outside the convex hull are null.

Maciek

P.S.
I'm very glad r.surf.nnbathy is usefull for you. Please drop me a line
if you publish anything where it's involved.

3

Maciej Sieczka wrote:

Luigi Ponti wrote:

  

I have successfully included natural neighbor interpolation in my 
mapping routine via r.surf.nnbathy: it seems slower (is that OK?)
than v.surf.idw but sounds a lot more robust from what I have read
around
    


I'd only like to stress that it really depends on what your input data
looks like and what the purpose is. Each interpolation method has it's
cons and pros. There is no single most robust method for all applications.

This is why I was worrying about the convex hull, please see
http://quartese.googlepages.com/idwvsnn

  

-- this is very useful if the interpolated raster is to be
used for further analysis, beyond visualization.

Currently, the interpolated raster output I get from r.surf.nnbathy
is confined to the convex hull delimited by the outermost input 
points/cells. Would there be a way to extend the interpolated raster
to the whole current region as in v.surf.idw
    


1. Note the hardcoded "-W 0" in the script. You can set it to something
smaller if you want to extrapolate beyond the convex hull. However, I
don't think this makes much sense - set "-W -2" and see how bad the
output DEM looks on the edges. And it can't look better, for how the
triangulation works. Triangulation is restricted to the convex hull
made of input points.

I have noticed that after sending the message: I tried with -1 and -3 and it really becomes weird. You are right that it does not make sense.

2. Even if your input raster is bigger than your current region,
r.stats, which is used in r.surf.nnbathy script to create the x,y,z
input for nnbathy, discards input raster's cells outside the current
region. Thus, the convex hull is created only for non-null input points
falling into the current region. And the convex hull shape will depend
on how the input non-null cells falling into the current region are
distributed. If you want a perfectly rectangular convex hull, make sure
there is a non-null input raster cell in each region's corner.

I'll think of automating this. Next vacation maybe :).

This would be great.

  

(in principle this
seems possible: outer Voronoi cells just get bigger)?
    


I don't exactly understand what you mean. Could you elaborate please?

What I meant is that the outer Voronoi cells may extend up to the region’s bounding box, with consequent new Voronoi cells for interpolation points outside the convex hull, please see
http://quartese.googlepages.com/voronoi

  

By looking at r.surf.nnbathy code and nnbathy --help, it seems that 
nnbathy should produce a rectangular grid of <cols> by <rows>. What
am I missing?
    


Yes, it does produce a rectangular grid. Only that all the cells
outside the convex hull are null.

Maciek

P.S.
I'm very glad r.surf.nnbathy is usefull for you. Please drop me a line
if you publish anything where it's involved.

Thanks Maciek for your answer and your work on r.surf.nnbathy.

4

Luigi Ponti wrote:

Maciej Sieczka wrote:

Luigi Ponti wrote:

(in principle this seems possible: outer Voronoi cells just get
bigger)?

I don't exactly understand what you mean. Could you elaborate
please?

What I meant is that the outer Voronoi cells may extend up to the
region's bounding box, with consequent new Voronoi cells for
interpolation points outside the convex hull, please see
http://quartese.googlepages.com/voronoi

nnbathy performs Delaunay triangulation, not Voronoi tesselation. I
(darkly) understand there is a mathematical relation between the two,
but I don't see how it could apply in your example. If this is just my
ignorance please explain more.

Anyway - what could be done to overcome the trimming to convex hull
effect: you need to make sure that the convex hull for your input cells
covers your whole area of interest - ie. enlarge region enough for
important input cells not to be omitted. If cells in the input raster
are distributed so that a convex hull covering the whole area of
interest is not possible, you need to enhance your input raster.

r.to.vect -z feature=point followed by v.hull might be faster for
testing than trial and error with r.surf.nnabthy itself.

Cheers,
Maciek

5

Maciej Sieczka wrote:

Luigi Ponti wrote:
  

Maciej Sieczka wrote:
    

nnbathy performs Delaunay triangulation, not Voronoi tesselation. I
(darkly) understand there is a mathematical relation between the two,
but I don't see how it could apply in your example. If this is just my
ignorance please explain more.

The Voronoi cells and Delaunay triangles are said to be ‘dual’ to one another and once one is known the other is completely defined (Sambridge et al. 1995). Natural neighbor interpolation relies on the construction of Voronoi diagrams for determining the areas used in computing weights. What Sakov did was to use the Triangle software to compute the underlying Delaunay triangulation (Fan et al. 2005). Or at least this is my understanding after a brief immersion into the available online refs! To be true, I haven’t been able to find out where Voronoi polygons show up in the nnbathy source code.

  • Sambridge M., Braun J., McQueen H., 1995. Geophysical parameterization and interpolation of irregular data using natural neighbours. Geophysical Journal International 122: 837-857. Available at http://wwwrses.anu.edu.au/geodynamics/nn/SBM95/SBM.html
  • Fan Q., Efrat A., Koltun V., Krishnan S., Venkatasubramanian S., 2005. Hardware-assisted natural neighbor interpolation. Workshop on Algorithm Engineering and Experiments (ALENEX). Available at www.cs.arizona.edu/~alon/papers/ast.ps
Anyway - what could be done to overcome the trimming to convex hull
effect: you need to make sure that the convex hull for your input cells
covers your whole area of interest - ie. enlarge region enough for
important input cells not to be omitted. If cells in the input raster
are distributed so that a convex hull covering the whole area of
interest is not possible, you need to enhance your input raster.

OK – got it: my input raster are weather stations, so it is going to be hard…

r.to.vect -z feature=point followed by v.hull might be faster for
testing than trial and error with r.surf.nnabthy itself.

Good advice.
Thanks again,
Luigi

6

Maciej Sieczka wrote:

Luigi Ponti wrote:
  

What I meant is that the outer Voronoi cells may extend up to the 
region's bounding box, with consequent new Voronoi cells for 
interpolation points outside the convex hull, please see 
[http://quartese.googlepages.com/voronoi](http://quartese.googlepages.com/voronoi)
    



Anyway - what could be done to overcome the trimming to convex hull
effect: you need to make sure that the convex hull for your input cells
covers your whole area of interest - ie. enlarge region enough for
important input cells not to be omitted. If cells in the input raster
are distributed so that a convex hull covering the whole area of
interest is not possible, you need to enhance your input raster.

Please, have also a look at
http://www.ems-i.com/gmshelp/Interpolation/Interpolation_Schemes/Natural_Neighbor_Interpolation.htm

This web page suggests an approach for extrapolation that would get a non-null value to all the points in the region: four additional input points are included in the NN interpolation that correspond to the corners of a user-defined region (a bounding box that includes all input points), by assigning them a value based on a different interpolation method (IDW). This way, the convex hull of input points has the same extent of the bounding box of the region. Also, I think that there is no need of giving the -W nnbathy flag a negative value and the resulting interpolated raster would still be inside the range of the input data because of the way IDW works.

It sounds feasible to add this functionality to r.surf.nnbathy, but I don’t know if the approach is correct from the GIS point of view (i.e., mixing two interpolation methods) or even more useful/better than just using IDW for the whole region.

Thanks and regards,

Luigi

7

Luigi Ponti wrote:

Maciej Sieczka wrote:

Luigi Ponti wrote:

What I meant is that the outer Voronoi cells may extend up to
the region's bounding box, with consequent new Voronoi cells for
interpolation points outside the convex hull, please see
http://quartese.googlepages.com/voronoi

Anyway - what could be done to overcome the trimming to convex
hull effect: you need to make sure that the convex hull for your
input cells covers your whole area of interest - ie. enlarge
region enough for important input cells not to be omitted. If
cells in the input raster are distributed so that a convex hull
covering the whole area of interest is not possible, you need to
enhance your input raster.

Please, have also a look at
http://www.ems-i.com/gmshelp/Interpolation/Interpolation_Schemes/Natural_Neighbor_Interpolation.htm

This web page suggests an approach for extrapolation that would get
a non-null value to all the points in the region: four additional
input points are included in the NN interpolation that correspond to
the corners of a user-defined region (a bounding box that includes
all input points), by assigning them a value based on a different
interpolation method (IDW). This way, the convex hull of input
points has the same extent of the bounding box of the region. Also,
I think that there is no need of giving the -W nnbathy flag a
negative value and the resulting interpolated raster would still be
inside the range of the input data because of the way IDW works.

It sounds feasible to add this functionality to r.surf.nnbathy, but
I don't know if the approach is correct from the GIS point of view
(i.e., mixing two interpolation methods) or even more useful/better
than just using IDW for the whole region.

It is possible to add such a functionality to r.surf.nnbathy, but I'd
rather leave it to the user to enhance his input raster the way he
prefers to, if necessary. Different interpolation methods or additional
data sources can be used to add the 4 corner points in question.

Cheers
Maciek

8

Maciej Sieczka wrote:

Luigi Ponti wrote:
  

Maciej Sieczka wrote:
    

Luigi Ponti wrote:
It sounds feasible to add this functionality to r.surf.nnbathy, but
I don't know if the approach is correct from the GIS point of view
(i.e., mixing two interpolation methods) or even more useful/better
than just using IDW for the whole region.
      




It is possible to add such a functionality to r.surf.nnbathy, but I'd
rather leave it to the user to enhance his input raster the way he
prefers to, if necessary. Different interpolation methods or additional
data sources can be used to add the 4 corner points in question.

Cheers
Maciek

Thanks Maciek,

Any pointer to getting the four corners? Maybe they are already there, i.e. the code that sets the working region for nnbathy?

set the working region for nnbathy (it’s cell-center oriented)

nn_n=echo $n | awk -v res="$nsres" '{printf "%.8f",$1-res/2.0}'
nn_s=echo $s | awk -v res="$nsres" '{printf "%.8f",$1+res/2.0}'
nn_w=echo $w | awk -v res="$ewres" '{printf "%.8f",$1+res/2.0}'
nn_e=echo $e | awk -v res="$ewres" '{printf "%.8f",$1-res/2.0}'

Though I don’t entirely understand what it does, it seems to be getting the coordinates of the center of cells located at the 4 corners of the region.

Thanks again and regards,

Luigi

9

Luigi Ponti wrote:

Maciej Sieczka wrote:

Luigi Ponti wrote:

Maciej Sieczka wrote:

Luigi Ponti wrote:

It sounds feasible to add this functionality to
r.surf.nnbathy, but I don't know if the approach is correct
from the GIS point of view (i.e., mixing two interpolation
methods) or even more useful/better than just using IDW for
the whole region.

It is possible to add such a functionality to r.surf.nnbathy, but
I'd rather leave it to the user to enhance his input raster the
way he prefers to, if necessary. Different interpolation methods
or additional data sources can be used to add the 4 corner points
in question.

Any pointer to getting the four corners? Maybe they are already
there, i.e. the code that sets the working region for nnbathy?

# set the working region for nnbathy (it's cell-center oriented)
nn_n=`echo $n | awk -v res="$nsres" '{printf "%.8f",$1-res/2.0}'`
nn_s=`echo $s | awk -v res="$nsres" '{printf "%.8f",$1+res/2.0}'`
nn_w=`echo $w | awk -v res="$ewres" '{printf "%.8f",$1+res/2.0}'`
nn_e=`echo $e | awk -v res="$ewres" '{printf "%.8f",$1-res/2.0}'`

Though I don't entirely understand what it does, it seems to be
getting the coordinates of the center of cells located at the 4
corners of the region.

It calculates the working region for nnbathy, so that it reflected the
current GRASS region. Since nnbathy is cell-center oriented each
N,S,E,W extent has to be shifted half a current GRASS region resolution
to avoid it's output missalignment with current GRASS region. This is
achieved with awk maths as above.

You can calculate the corners of your current GRASS region from the
output of g.region -g.

Maciek

10

.

Since nnbathy is cell-center oriented each
N,S,E,W extent has to be shifted half a current GRASS region resolution
to avoid it's output missalignment with current GRASS region.

I always though that GRASS was cell-center oriented..

Carlos

--
+-----------------------------------------------------------+
              Carlos Henrique Grohmann - Guano
  Visiting Researcher at Kingston University London - UK
  Geologist M.Sc - Doctorate Student at IGc-USP - Brazil
Linux User #89721 - carlos dot grohmann at gmail dot com
+-----------------------------------------------------------+
_________________
"Good morning, doctors. I have taken the liberty of removing Windows
95 from my hard drive."
--The winning entry in a "What were HAL's first words" contest judged
by 2001: A SPACE ODYSSEY creator Arthur C. Clarke

Can't stop the signal.

11

Maciek wrote

Carlos "Guâno" Grohmann wrote:

Since nnbathy is cell-center oriented each
N,S,E,W extent has to be shifted half a current GRASS region resolution
to avoid it's output missalignment with current GRASS region.

I always though that GRASS was cell-center oriented..

What do you mean?

Explanation of why the shift in r.surf.nnbathy in necessary:

The region:

$ g.region -g
n=404456
s=404453
w=315942
e=315945
nsres=1
ewres=1
rows=3
cols=3
cells=9

The input raster (* stands for null):

315942.5 404455.5 20
315943.5 404455.5 74
315944.5 404455.5 43
315942.5 404454.5 29
315943.5 404454.5 *
315944.5 404454.5 65
315942.5 404453.5 25
315943.5 404453.5 25
315944.5 404453.5 49

If the half-a-resolution shift IS applied (nnbathy -W 0 -P alg=nn -n
3x3 -x 315942.5 315944.5 -y 404455.5 404453.5), as described in
previous post, the output of r.surf.nnbathy in=xxx out=yyy is alligned
to the input raster:

315942.5 404455.5 20
315943.5 404455.5 74
315944.5 404455.5 43
315942.5 404454.5 29
315943.5 404454.5 48.25
315944.5 404454.5 65
315942.5 404453.5 25
315943.5 404453.5 25
315944.5 404453.5 49

If the shift IS NOT applied (nnbathy -W 0 -P alg=nn -n 3x3 -x 315942
315945 -y 404456 404453), the output is missaligned and interpolated
wrong, due to nnbathy working region set wrong:

315942 404456 NaN
315943.5 404456 NaN
315945 404456 NaN
315942 404454.5 NaN
315943.5 404454.5 48.25
315945 404454.5 NaN
315942 404453 NaN
315943.5 404453 NaN
315945 404453 NaN

Maciek

12

OK. thanks for the expalnation.

carlos

On 9/6/07, Maciej Sieczka <tutey@o2.pl> wrote:

>> Maciek wrote
> Carlos "Guâno" Grohmann wrote:

>> Since nnbathy is cell-center oriented each
>> N,S,E,W extent has to be shifted half a current GRASS region resolution
>> to avoid it's output missalignment with current GRASS region.

> I always though that GRASS was cell-center oriented..

What do you mean?

Explanation of why the shift in r.surf.nnbathy in necessary:

The region:

$ g.region -g
n=404456
s=404453
w=315942
e=315945
nsres=1
ewres=1
rows=3
cols=3
cells=9

The input raster (* stands for null):

315942.5 404455.5 20
315943.5 404455.5 74
315944.5 404455.5 43
315942.5 404454.5 29
315943.5 404454.5 *
315944.5 404454.5 65
315942.5 404453.5 25
315943.5 404453.5 25
315944.5 404453.5 49

If the half-a-resolution shift IS applied (nnbathy -W 0 -P alg=nn -n
3x3 -x 315942.5 315944.5 -y 404455.5 404453.5), as described in
previous post, the output of r.surf.nnbathy in=xxx out=yyy is alligned
to the input raster:

315942.5 404455.5 20
315943.5 404455.5 74
315944.5 404455.5 43
315942.5 404454.5 29
315943.5 404454.5 48.25
315944.5 404454.5 65
315942.5 404453.5 25
315943.5 404453.5 25
315944.5 404453.5 49

If the shift IS NOT applied (nnbathy -W 0 -P alg=nn -n 3x3 -x 315942
315945 -y 404456 404453), the output is missaligned and interpolated
wrong, due to nnbathy working region set wrong:

315942 404456 NaN
315943.5 404456 NaN
315945 404456 NaN
315942 404454.5 NaN
315943.5 404454.5 48.25
315945 404454.5 NaN
315942 404453 NaN
315943.5 404453 NaN
315945 404453 NaN

Maciek

--
+-----------------------------------------------------------+
              Carlos Henrique Grohmann - Guano
  Visiting Researcher at Kingston University London - UK
  Geologist M.Sc - Doctorate Student at IGc-USP - Brazil
Linux User #89721 - carlos dot grohmann at gmail dot com
+-----------------------------------------------------------+
_________________
"Good morning, doctors. I have taken the liberty of removing Windows
95 from my hard drive."
--The winning entry in a "What were HAL's first words" contest judged
by 2001: A SPACE ODYSSEY creator Arthur C. Clarke

Can't stop the signal.

13

Carlos wrote:

> Since nnbathy is cell-center oriented each
> N,S,E,W extent has to be shifted half a current GRASS region
> resolution to avoid it's output missalignment with current GRASS
> region.

I always though that GRASS was cell-center oriented..

the region is set from the outer bound of a cell.

For example, a 1x1 cube is centered at 0.5,0.5 within a region of 0,1 0,1:

G63> g.region -p n=1 s=0 e=1 w=0 res=1
G63> r.mapcalc one=1
G63> r.out.xyz one
0.5|0.5|1

see the first few lines of rasterintro.html

Hamish