Hi Rupert,

I'm replying to your two last posts.

I found these two lines in the InterpolateDEM function:

ix = INT((x-xb0)/dbx)+1

iy = INT((y-yb0)/dbx)+1

Surely the second line should be divide by dby not dbx? This would be a problem for grids with non-square cells with variation in the y-direction.

Yes, this is indeed a bug. Thanks for seeing it!

The default in the 2d grid interpolator in the event of one to three neighbouring points having the noData value is to take the mean of valid points. I have written a few lines to instead carry out linear interpolation. Shall I set this as the default (which makes more sense to me) or as an option (for full backwards compatibility)?

I'm not sure a linear interpolation is possible in all cases:

- with only one node, it can only be estimated as the value at that node

- with two nodes, it should work, but one as to determine the direction of the interpolation (x or y)

- with three nodes, it would work only if the point at which you want the value is included in the triangle defined by the three nodes.

So, I would be interested to know more about the strategy you have developed to handle this interpolation.

An other solution should be to use extrapolation from the neighbouring cells.

Nevertheless, missing data should only concern a very limited number of points and at the margin of the domain. Else, you are in trouble...

Regards

Olivier