[GRASS-dev] new flag to interpolate from nearby raster cells for v.what.rast

GRASS GIS GRASS GIS developer
1

Hi,

I just added in the dev branches two new flags for v.what.rast. The first
is a flag to print the query results to stdout instead of updating a
table, since for a large vector points map the table update is
effectively unusable unless you let it run overnight (bug #2131).

The second flag changes the default containing-grid-cell method to a
weighted avg. of the 4 nearest raster cell centers (IDW). The search radius
is not very big, which does lead to some cell edge artifacts if you look
closely, but it runs almost at the same speed as the nearest-neighbour
query and is not overly complicated, so I'd hope "good enough". See the
attached screenshot and try to spot the flaws (becomes obvious in NVIZ).

This is mostly useful when going from a very coarse resolution raster map
to a very dense set of vector points and you want to avoid large steps in
value at the cell edges. A common occurrence of this is going from world
or continental scale climatic or weather raster data -> input to a local
finite element mesh-grid*, and the big jumps in value risk causing
shockwaves or other unpleasant effects in the physics of the finite
element model.

[*] e.g. see "Gmsh" http://geuz.org/gmsh/gallery/MURUROA.png

The same technique could be applied to r.what if there was the need, Radim
adapted the v.what.rast module from it so the code bases are pretty similar.
The method seems to work pretty well, but I'd welcome a better coder than
myself to continue improving it. I considered parallelizing it with OpenMP,
but it will run through a million points in ~4 seconds, so perhaps not
worth the trouble. Another idea is to allow it to work without topology,
so fast raster querying from massive point swarms (e.g. lidar).

please test!

enjoy,
Hamish

(attachments)

v.what.rast_interp.png

2

On Thu, Nov 14, 2013 at 11:15 AM, Hamish <hamish_b@yahoo.com> wrote:

Hi,

I just added in the dev branches two new flags for v.what.rast.
[...]

The second flag changes the default containing-grid-cell method to a
weighted avg. of the 4 nearest raster cell centers (IDW). The search radius
is not very big [...]

AFAIK, IDW takes as argument the radius which is fixed to 1 which is
questionable. Since v.what.rast samples a raster value at a given
location, why not use the standard raster sampling methods
                nearest: nearest neighbor
                linear: linear interpolation
                cubic: cubic convolution
                lanczos: lanczos filter
                linear_f: linear interpolation with fallback
                cubic_f: cubic convolution with fallback
                lanczos_f: lanczos filter with fallback

which are used by r.proj and i.rectify? The concept is essentially the
same: for a given location, estimate the corresponding raster value.

Markus M

3

Hi,

Hamish wrote:

I just added in the dev branches two new flags for v.what.rast.
[...]

The second flag changes the default containing-grid-cell method to a
weighted avg. of the 4 nearest raster cell centers (IDW). The search radius
is not very big [...]

AFAIK, IDW takes as argument the radius which is fixed to 1 which is
questionable.

I'm not sure I follow what you mean, but the "distance" in the 4-way inverse
distance weighting in the new v.what.rast -i interpolate flag is calculated
for each of the closest four raster cells' center coord by running
G_distance() vs. the vector point's position, and then the weighting in the
weighted average for each of those four cells' data values is by 1/distance^2.

It doesn't assume that the 4 surrounding cells are equally distant:

weightsum = valweight = 0;

weight[i] = 1.0 / (distance[i] * distance[i]);
weightsum += weight[i];

valweight += weight[i] * nearby_d_val[i];
value = valweight / weightsum; // 0./0. == nan

The other thing of note is that if any of the nearest four cells are null
then it discards it from the weighted average, and proceeds with n-1 elements.
I find this a preferable compromise between the null if-any-nulls and the
fallback to nearest neighbor if-any-nulls methods below; YMMV.

Since v.what.rast samples a raster value at a given
location, why not use the standard raster sampling methods
nearest: nearest neighbor
linear: linear interpolation
cubic: cubic convolution
lanczos: lanczos filter
linear_f: linear interpolation with fallback
cubic_f: cubic convolution with fallback
lanczos_f: lanczos filter with fallback

which are used by r.proj and i.rectify? The concept is essentially the
same: for a given location, estimate the corresponding raster value.

what I've called IDW is similar to the bilinear method above, with the
exception of how surrounding NULLs are handled, and nearest is similar
to the default behavior of v.what.rast. If someone wants to take it up
and expand the options to use Rast_interp_bicubic() & co., I'm happy for
it, but would like to preserve the null handling and speed.

The main bug so far is if the print-only flag is used, it doesn't yet
bypass the checks for if the vector map has a table, is in the current
mapset (is alterable), and to remove the need for topology.

regards,
Hamish