1

I have a set of GPS stations, whose coordinates I know (x,y,z), and for each station, I also have an error (e). These stations are of course unevenly spaced, otherwise it would be too easy. The thing is, in order to compute my error e for a station, I have only used said station, but I also want to take the others into account.

My problem is: given a set of unevenly spaced (x,y,z,e) points, how can I interpolate e in function of the spatial distance between points? The interpolation would need not to be exact, since I am recalculating e on points where I already have it. Also, I am looking for something cleaner than inverse distance or likewise stuff. Splines would be nice, for example.

From what I have read, the splev function of the scipy.interpolate package seems to do the trick, but I cannot understand how it works or what I am supposed to give it as arguments.

I someone could either explain how this is supposed to work, or point me to another method, it would be great.

Bach
  • 6,145
  • 7
  • 36
  • 61
Efferalgan
  • 1,681
  • 1
  • 14
  • 24
  • Why do you think that "inverse distance or likewise stuff" isn't a "clean" solution? Also, as DrV asks below, do you truly have 3D points, or can you get away with 2.5D (Earth's surface plus elevation)? – Daniel Pryden Jun 15 '14 at 22:35
  • I need to really smoothe my errors, and no to get a too "spiky" solution. They really are correlated, and there is some noise on top of that. – Efferalgan Jun 15 '14 at 23:16

1 Answers1

1

If I understood your question correctly, you have a point x,y,z in space, and you would like to calculate the error by interpolating from the known stations. You also suggest that the validity of an error reading depends on the distance from the point where the error is known.

So, for a point x,y,z you can calculate its distance from each of the known receiving stations. Then you have some weight function calculated from each of these distances. Finally, you take a weighted average by the weight functions (or possibly do some other trickery to eliminate outliers).

How about this:

# known station coordinates (number of rows matching number of stations)
coords = array([
    (x1, y1, z1),
    (x2, y2, z2),
    ... ])
# respective error values (number of items matching number of stations)
err_values = array([
    e1,
    e2),
    ... ])

# p is a three-element array representing our coordinates
# distances will contain the spatial euclidian distances to the stations 
distances = numpy.linalg.norm(coords - p[None,:], axis=1)

# get the weights somehow
weights = my_weight_function(distances)

# return the weighted average
return numpy.average(err_values, weights=weights)

There is one more trick which might be useful especially in this case. The last statement could be replaced by:

return numpy.sum(err_values * weights) / (eps + numpy.sum(weights))

Basically a weighted sum, but a small number eps added to the denominator. The point here is that as we talk about an error, it should be zero very far away from the known points. Otherwise we would often have an average of the known errors as the error on the other side of the globe, which is not reasonable. The only reasonable assumption is that the error is zero far away. It is not, but we do not know any better, and thus zero is the best guess.

If I understood your problem in a wrong way, let me know. (If you are thinking of the interpolation problem as a way to increase accuracy on the surface of the earth, you actually have a 2d problem on the surface of the globe, not a real 3D problem.)

DrV
  • 22,637
  • 7
  • 60
  • 72
  • Isn't this just the IDW (inverse-distance-weighted) algorithm? – Daniel Pryden Jun 15 '14 at 22:37
  • Well, it is the *DW algorithm, you choose the weights :) I suspect that trying to find anything more intelligent would require much more knowledge on the error source. Just using splines does not help much, if the underlying problem does not produce splines. The distance-weighted algorithms are at least reasonably stable. – DrV Jun 15 '14 at 22:43
  • Thank you, but that is basically IDW, and I cannot use that. There is a strong correlation between my errors, and I suspect the underlying problem is more likely to produce spline-like shapes than things with spikes like what IDW usually produces. Moreover, like I said, my stations are unevenly spaced, with near-empty areas and clusters, which will attract the corrections toward them. I could maybe get away with Earth surface + elevation, but that would require a huge transformation of my data (which is quite important), which is in metres from the centre of the Earth. – Efferalgan Jun 15 '14 at 23:07
  • The idea here is that you do not use IDW (as it is really a spiky beast), but you use something else than 1/x as the weight function. Your weight function should reflect the correlation between stations as a function of distance from the station. I would experiment with, e.g., a gaussian with different widths to see if it is useful. The problem with, e.g., 3D splines is that your data is locally 2D data ("Earth is flat"), and fitting something 3D to that data easily produces very odd results in the third dimension. – DrV Jun 16 '14 at 05:54
  • My stations are scattered across all New-Zealand, so I have managed to bypass the problem by finding the right projection from the NZ surveying institute, and I am eventually working in 2D. But just as a theoretical problem, I was wondering how to do 3D splines with Python. – Efferalgan Jun 16 '14 at 23:05
  • Thank you for accepting my answer despite it being exactly what you did not want to have! As your data is very flat (I know NZ is not that flat), you cannot really use any 3D interpolation, as they would become extrapolations. I have actually never used 3D splines, as I have never had enough evenly-spread data points to create a reasonable fit. Interpolations are only valid when they are "in" the set of known data points. – DrV Jun 17 '14 at 09:03