Geographic Interpolate in Python

Question

I have 500 points with longitude x, latitude y, altitude z, and the value at these points.

On the other hand, I have other 200 points than I would like to interpolate, where the latitude, longitude, and altitude are known.

I would like to interpolate considering the altitude of the points and geography between these points, my maps are Spain.

An example:

I have 2 points (x,y,z) and (x',y',z'). The real distance is black line, a polynomial interpolation is blue line (approx. distance), and euclidean distance between these 2 point is red line. I would like to obtain the blue or black line distance.

The following example which takes (x,y) as input would fit:

https://nbviewer.jupyter.org/github/cjohnson318/geostatsmodels/blob/master/notebooks/KrigingExample.ipynb

but I would like to also manage the altitude z as a input parameter.

Some libraries in Python? Some tutorial?

Answer 1

I'll try to answer your question using OpenTURNS platform.

Let's consider that Spain is a square 1000 x 1000 km and that your 500 points are randomly distributed over the surface

import openturns as ot
import numpy as np

# initiate a sample of size 500 with 2 coordinates
inputdata = ot.Sample(500, 2)
# 1st column random between 0 and 1000
inputdata[:,0] = ot.Uniform(0,1000).getSample(500)
# 2nd column random between 0 and 1000  
inputdata[:,1] = ot.Uniform(0,1000).getSample(500)

Then let's assign a height for each of these points. OpenTURNS allows to define symbolic functions :

height = ot.SymbolicFunction(["x","y"], ["10 +10 * (x + y) / 1000 + 10 * ((x + y) / 1000) * sin( 3 * x * pi_ / 1000 )*cos(5 * y * pi_ / 1000)"])
outputdata = height(inputdata)

Now we would like to interpolate the data to estimate the height of any point on the map. The Kriging method allows to do so but you'd better know some information about your problem (general trend, correlation between the heights of 2 distant points).

# dimension of the input data
dimension = 2
basis = ot.ConstantBasisFactory(dimension).build()
covarianceModel = ot.SquaredExponential(dimension)

Then we just call the kriging algorithm to do the interpolation

algo = ot.KrigingAlgorithm(inputdata, outputdata, covarianceModel, basis)
algo.run()
result = algo.getResult()
metamodel = result.getMetaModel()

metamodel is exactly the function you want!

# gives the inferred height of the point (x = 123, y = 967)
metamodel([123, 967])
>>> [12.2225]

Should you like to draw the result, you can calculate the predicted values on a meshgrid of your square

gridx = np.arange(0.0,1001,10)
nx = len(gridx)
gridy = np.arange(0.0,1001,10)
ny = len(gridx)
X, Y = np.meshgrid(gridx, gridy)

predictions = np.array(metamodel([[xi,yi] for (xi, yi) in zip(X.ravel(),Y.ravel())])).reshape(nx,ny)

then you can use matplotlib to view the result:

import matplotlib.pylab as plt
plt.figure()
vmin = predictions.min()
vmax = predictions.max()
plt.pcolor(X, Y, predictions, cmap='viridis', vmin=vmin, vmax=vmax)
plt.scatter([d[0] for d in inputdata], [d[1] for d in inputdata],  c = [d for d in outputdata], s=2, edgecolor = "white", cmap='viridis', vmin=vmin, vmax=vmax)
plt.colorbar()
plt.show()

You can also view it in 3D :-)

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
    
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, predictions, cmap=cm.viridis,
                           linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)
    
plt.show()