[GRASS-user] Raster Hole/NODATA filling programs

GRASS GRASS user
1

Folks,

I have a client interested in an alternative implementation of:

   http://www.pcigeomatics.com/cgi-bin/pcihlp/GRDINT

This is a program for filling nodata areas in a raster by interpolating
in from the edges. The GRDINT program itself is primarily targetted at
interpolating from rasterized contour data, but it is also useful in a wide
variety of other contexts - to fill missing data holes in raster data.

I've looked at r.fillnulls and r.surf.idw and they both seems to attack
aspects of this problem. In the case of r.surf.idw it talks about
interpolating from the "n nearest points" but this seems entirely unsuitable
for filling in large nodata areas since it will tend to have a substantial
discontinuity near the middle of the filled region as the n-nearest points
flip from top to bottom or left to right.

The r.fillnulls *appears* to scan the entire boundary of a nodata area to
build a point list of boundary values used for spline based interpolation.
I'm concerned this would not work well for situations like contour data with
large nodata areas or for situations where relatively sparse data has
been applied to a raster and we need to interpolate from that (so essentially
most of the image might be one connected nodata area!).

Are there other GRASS algorithms that I should be looking at for this task?

Am I misconstruing the weaknesses of r.fillnulls and r.surf.idw? Perhaps
it is somewhat foolhardy for me to want a one-size-fits-all nodata filler?

If I implemented an algorithm like GRDINT (possibly lacking the morphological
aspects, but based on 8-direction searches) would it be a useful addition to
GRASS?

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | President OSGeo, http://osgeo.org

2

On May 14, 2008, at 2:25 PM, Frank Warmerdam wrote:

The r.fillnulls *appears* to scan the entire boundary of a nodata area to
build a point list of boundary values used for spline based interpolation.
I'm concerned this would not work well for situations like contour data with
large nodata areas or for situations where relatively sparse data has
been applied to a raster and we need to interpolate from that (so essentially
most of the image might be one connected nodata area!).

Yeah, r.fillnulls has problems with large holes. I used it on the SRTM data, and large holes get block anomalies, probably because of the SEGMAX (for v.surf.rst) used by r.fillnulls. I played with that a bit in a customized r.fillnulls, but the bigger the SEGMAX the longer the processing, and it only helped a little.

I didn't pursue any further improvements to fillnulls, though I did experiment with filling at a coarser resolution, converting to points and using v.surf.rst with those as an enhancement to a similar process to r.fillnulls (it gave the full resolution v.surf.rst a little more interior meat to work from). But I got distracted by "real" work and dropped the experiment.

-----
William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>
http://www.kyngchaos.com/

[Trillian] What are you supposed to do WITH a maniacally depressed robot?

[Marvin] You think you have problems? What are you supposed to do if you ARE a maniacally depressed robot? No, don't try and answer, I'm 50,000 times more intelligent than you and even I don't know the answer...

- HitchHiker's Guide to the Galaxy

3

Frank:

> The r.fillnulls *appears* to scan the entire boundary
> of a nodata area to build a point list of boundary
> values used for spline based interpolation.

yes, it is good for filling small holes.

> I'm concerned this would not work well for situations like contour
> data with large nodata areas or for situations where relatively
> sparse data has been applied to a raster and we need to interpolate
> from that (so essentially most of the image might be one connected
> nodata area!).

probably not. In the case of contour lines have a look at the r.surf.contour module, it seems to do a nice job there. In the case of sparse data starting points, v.surf.rst will do a nice job. The method to use depends heavily on the input data, I'm not sure if you will find a one size fits all solution..

William:

Yeah, r.fillnulls has problems with large holes. I used it
on the SRTM data, and large holes get block anomalies, probably
because of the SEGMAX (for v.surf.rst) used by r.fillnulls.
I played with that a bit in a customized r.fillnulls, but the bigger
the SEGMAX the longer the processing, and it only helped a little.

I didn't pursue any further improvements to fillnulls, though I
did experiment with filling at a coarser resolution, converting
to points and using v.surf.rst with those as an enhancement to a
similar process to r.fillnulls (it gave the full resolution
v.surf.rst a little more interior meat to work from).

Right. v.surf.rst has problems when the density of the data points changes suddenly. (such as a dense ring of points with a huge gap in the middle)
If you have randomly spaced starting points of a constant density, v.surf.rst (tps) is a hard thing beat.
see http://grass.osgeo.org/wiki/RST_Spline_Surfaces#Troubleshooting

I recently updated the v.random.cover module in the wiki addons for that task William, ie create a idealized surface from just the loose density points using v.surf.rst + use v.hull (or if using r.surf.nnbathy it does the hull cropping automatically), then use v.random.cover to create n random points and take the elevation from the newly interpolated helper map. Patch the supplimental pseudo-points into your original dataset and then run v.surf.rst without the problems cause by the starting points density change.

It would be interesting to explore radial basis functions like r^2log(r) or r^4log(r) with continous derivates to get away from the poorly patched segment problems. But maybe that doesn't scale well past 10k start points?
or maybe that is the same as segmax=n starting points?

other modules which may be of interest,

* r.surf.contour: sounds like what you are looking, but be careful that it needs an update for floating point output and not to treat 0 as no-data.
To work around this in the past I have multiplied the elevation map by 100000 before starting, then reclassifyling 0 as 1. After processing I divide the resulting map by 100000. The search time can be long if the starting image includes big gaps, but the results are rather nice.

* see also r.surf.idw2. Be careful with these modules, they may be on the "waiting for FP update" queue as well; I'm not sure.

* r.surf.nnbathy from the wiki addons:
  http://grass.osgeo.org/wiki/GRASS_AddOns#r.surf.nnbathy
(nn= natural neighbor) It does a nice job and offers a few different methods. It is well worth the time spent to have a look at.
I'm not sure about its FP/int status.

Hamish