search

harphub / meteogrid

R package for working with gridded meteorological data
https://harphub.github.io/meteogrid/
Other
0 stars 6 forks source link

Bug in point.interp #11

Open FlorianW-ZAMG opened 3 years ago

FlorianW-ZAMG commented 3 years ago

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 _transformgeofield.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 commented 3 years ago

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.