Bug in point.interp

[FlorianW-ZAMG]

I assume there is a bug in point.bilin function in interpol.R. When I do gribfile extraction with harp using bilinear interpolation this routine is called from transform_geofield.R. The point.bilin function returns a n x n matrix with n equals number of stations. The result variable in point.bilin looks like:

[1] "RESULT"
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 289.2113 288.8490 286.9476 288.1653 288.5911 287.0800 287.5619 289.2113
[2,] 287.3979 286.7967 287.1780 285.0586 284.4355 288.6827 288.4766 287.3979
[3,] 288.0519 287.9129 287.4847 286.3963 287.5759 287.6278 288.8027 288.0519
[4,] 287.6618 287.0082 286.4630 284.6508 287.5166 287.0171 287.5300 287.6618
....

I do not know much about harp but binding the indices like:
result <- weights$w00*infield[cbind(weights$i0, weights$j0)] + weights$w01*infield[cbind(weights$i0, weights$j1)] + weights$w10*infield[cbind(weights$i1, weights$j0)] + weights$w11*infield[cbind(weights$i1, weights$j1)]

seems to fix this. The result now looks like:

[1] "RESULT"
[1] 289.2113 286.7967 287.4847 284.6508 284.2716 287.8197 288.6032 286.6491
[9] 284.7303 285.9676 285.4736 290.3426 288.6061 285.8485 288.5417 290.1415
[17] 289.0653 292.6197 292.1076 292.1517 289.9855 291.0373 290.3357 290.4631
...

However there might be a more 'R-like' solution to this issue, so I did not commit a pull request.

Best regards,
Florian

[adeckmyn]

You are right. I'm afraid I introduced this bug in the latest update. I actually fixed another (worse) bug in the interpolation, but I guess I only checked my clean-up with single points and full 2d fields, not with a list of stations (n>1). The older version actually used cbind more or less in the way you propose. I'll update the package asap. Thanks for reporting this.