Interpolation of non-gridded data in 2 and 3 dimensions

Announcements
1

I have published a very preliminary project on Github - GitHub - arjenmarkus/interpolation2d3d: Interpolation in 2D and 3D, object-oriented interfaces - that modernises several classical packages available on Netlib for interpolating spatial data. As @ivanpribec pointed out, I need to clarify the actual licensing conditions, as it uses code from ACM TOMS. I am working on that :slight_smile: Meanwhile, have a look at the interface I propose - it hides much of the complexity of the original code, thanks to the object-oriented features. (It does lack more elaborate examples and tests, so that is to be considered ā€œwork-in-progressā€.)

2 Likes
2

@Arjen, you might also have a look at Alexander Barthā€™s ā€œoptimalā€ interpolation code (based on Kriging) in the NASA CFDtools github repository.

From the README

ā€˜ā€˜This Optimal Interpolation module provides a powerful means of treating
unstructured datasets in n-space, whether they be clouds of points from
laser scans of some objects, irregularly spaced function values that need
interpolating to new target independent variables, or whatever. It belongs
in the family of methods best known as Kriging, which determine the linear
combination of neighboring data points that minimizes each interpolation
error in some sense under certain statistical assumptions (hence the use
of the term optimal).ā€™ā€™

The license is BSD and NASAā€™s open source license

1 Like
3

Thanks for this reference. I did not know it. It does seem of a completely different class of algorithms, but it should be useful as part of an interpolation toolset.

4

My understanding of Kriging (and I could be wrong) based interpolation methods is that they are in the same family as Radial Basis Functions (RBFs). If you have time you might also take a look around for some RBF packages and associated literature.

5

Had only a brief look at the tests, but would it make sense for the ā€œinitializerā€ to be more of a functional style? I.e. instead of

call shepard%create_mesh( x, y ,f )

do

shepard = qshep_2d(x, y, f)

Also, and I havenā€™t run it so Iā€™m not sure what it does now, but what is returned when an error occurs? I.e.

    px = 1.5
    py = 1.5

    write(*,*) px, py, shepard%interpolate_point( px, py )
    call shepard%error_status( error_code, error_message )
    write(*,*) 'Error: ', error_code, ' -- ', trim(error_message)

seems to suggest that an error is expected, but what would prompt a new user to look for and find the error_status procedure? What happens if I perform multiple interpolations in an expression? Like

foo = shepard%interpolate_point(1.25, 1.25) + shepard%interpolate_point(0.75, 0.75)

The standard says nothing about which is evaluated first, and so am I guaranteed to receive an error if one of them fails? What error do I get if both fail?

I commend the effort. This is the kind of library that Iā€™m sure will prove to be quite useful in many applications.

6

Thanks for your comments - I have thought of an initialisation function, but in the end chose this style. With a function you would have to copy the data on assignment, so the style I chose might be more efficient. Alternatively, you can define the appropriate assignment procedure ā€¦

I will look into this after my holidays :slight_smile:

7

everythingfunctional has given us one good reason for using a subroutine instead of a function: if a statement like foo = f(x1) + f(x2) fails because of something done inside the function f, it might be hard to say whether it was f(x1) or f(x2) that caused the trouble. Another good reason that hit me yesterday was that if print *,f(x) did not fail but there was another print statement inside f, the output looked somewhat garbled. So I turned my function f into a subroutine with one intent(out) argument.

8

In the last example, where f(x) is in an i/o list and itself prints output, that case is specifically not allowed by fortran. Up until recently, no i/o of any kind at all could be done in this manner. More recently, some exceptions are allowed, such as doing i/o to separate output units.

9

I know that if a function is in an output list to unit u and itself writes to the same unit u then that is forbidden. However both gfortran and ifort let me do it with the program below. That is not a bug in either compiler because it is one of the many errors that the standard does not require to be diagnosed. It was revealed (as often in such cases) by the different compilers giving different outputs.

program badwrite ! Forbidden by the Fortran standard
  implicit none
  write(*,"(A,I0,A) ") 'f()="',f(),'"'
contains
  integer function f()
    f = 666
    write(*,"(A)") 'Hello world!'
  end function f
end program badwrite