Interpolate undulations for given points using a geoid model

Question

I have two databases, the first one (F.tsj) contains gps points: coordinates and ellipsoid heights and the second one is simply a geoid model where every line stores latitude, longitude, and undulation. My simple program retrieves four adjacent points from the geoid grid database for every point stored in the first databse. Both databases store coordinates in decimal degrees.
Here's my code to build geoid file database from csv file:

import csv,sqlite3
'''this script inserts data stored in csv file into sqlite database'''

coordinates= csv.reader(open("geoid.csv"))

database = sqlite3.connect("geoid.db")
cursor = database.cursor()
database.execute("create table coords(name, lat INTEGER, lon INTEGER, und INTEGER);")
database.executemany("insert into coords(name,lat,lon,und) values (?, ?, ?, ?)",coordinates)
database.commit()
database.close()

Here's example input csv file, that contains ID,Latitude,Longitude,Undulation:
000001,54.9666669,18.7,28.442
000002,54.9166669,18.6833331,28.512
000003,54.9,18.6833331,28.523
000004,54.8833331,18.6666669,28.572
000005,54.8666669,18.6666669,28.585
000006,54.85,18.6666669,28.6
000007,54.8333331,18.0666669,29.891
000008,54.8333331,18.0833331,29.855
000009,54.8333331,18.1833331,29.637
000010,54.8333331,18.25,29.488
000011,54.8333331,18.2666669,29.45
000012,54.8333331,18.2833331,29.413
000013,54.8333331,18.3,29.376
This is numeric geoid model stored in ASCII csv file.

Now I would like to find four adjacent points for our example point:

import sqlite3
#find four adjacent points
def find_adjacent_coords(db, lat, lon, step=0.041667):
    coords_range = lat-step, lat+step, lon-step, lon+step
    return db.execute("""select lat, lon, und from coords where
lat > ? and lat < ? and
lon > ? and lon < ?""", coords_range).fetchall()
#connect database with gps points
database = sqlite3.connect('F.tsj')
cursor = database.cursor()
cursor.execute("select C1, C2, C3 from tblSoPoints")
#C1 is latitude
#C2 is longitude
#C3 is ellipsoid height
results = cursor.fetchall()
for line in results:
    B = line[0]
    L = line[1]
    ellps_height = line[2]
#connect to geoid database
    db = sqlite3.connect('geoid.db')
    n = (find_adjacent_coords(db, a, b))

Example gps point:

54.4786674627, 17.0470721369, 86.5003132338

and adjacent points from geoid:

[(54.5, 17.041667, 31.993), (54.5, 17.083333, 31.911), (54.458333, 17.041667, 31.945), (54.458333, 17.083333, 31.866)]<br>

Then I would like to interpolate undulation for our point using these four adjacent points with known undulations to calculate orthometric height.

Edit

Solution No 1

from __future__ import division

def bilinear_interpolation(x, y, points):
    points = sorted(points)               # order points by x, then by y
    (x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points

    if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2:
        raise ValueError('points do not form a square')
    if not x1 <= x <= x2 or not y1 <= y <= y2:
        raise ValueError('(x, y) not within the square')

    return (q11 * (x2 - x) * (y2 - y) +
            q21 * (x - x1) * (y2 - y) +
            q12 * (x2 - x) * (y - y1) +
            q22 * (x - x1) * (y - y1)
           ) / (x2 - x1) * (y2 - y1)

n = [(54.5, 17.041667, 31.993),
         (54.5, 17.083333, 31.911),
         (54.458333, 17.041667, 31.945),
         (54.458333, 17.083333, 31.866),
]

B = 54.4786674627
L = 17.0470721369
print bilinear_interpolation(B, L, n)
import doctest
print doctest.testmod()

...and here's the solution No 2:

from __future__ import division
import doctest
def bilinear_interpolation(x, y, points):

# http://en.wikipedia.org/wiki/Bilinear_interpolation
    points.sort()               # order by x, then by y
    xs, ys, zs = zip(*points)
# verify it is a square
    assert xs[0] == xs[1] and xs[2] == xs[3]
    assert ys[0] == ys[2] and ys[1] == ys[3]
    x1 = xs[0]
    x2 = xs[2]
    y1 = ys[0]
    y2 = ys[1]
# verify interpolation versus extrapolation
    assert x1 <= x <= x2, (x1, x, x2)
    assert y1 <= y <= y2

    q11, q12, q21, q22 = zs
    return (q11 / ((x2 - x1) * (y2 - y1)) * (x2 - x) * (y2 - y) +
            q21 / ((x2 - x1) * (y2 - y1)) * (x - x1) * (y2 - y) +
            q12 / ((x2 - x1) * (y2 - y1)) * (x2 - x) * (y - y1) +
            q22 / ((x2 - x1) * (y2 - y1)) * (x - x1) * (y - y1)
           )
B = 54.5
L = 17.061667

n = [(54.5, 17.041667, 31.993),
     (54.5, 17.083333, 31.911),
     (54.458333, 17.041667, 31.945),
     (54.458333, 17.083333, 31.866),
]
print doctest.testmod()

print bilinear_interpolation(B, L, n)

Credits for these solutions goes to:

Raymond Hettinger

Answer 1

You have two tasks: (a) find the grid data for the vertices of a grid cell in which a GPS location (x,y) lies and (b) interpolate those data.

An efficient way to accomplish (a) is to store the grid as an array with a single index. Often grids are stored by rows with data progressing left to right within each row. A grid is determined by the coordinates of its origin, (x0,y0) (which will be (-180, -90) for a worldwide grid in decimal degrees), its cellsize c (which is 1/24 in your case), and the number of points per row (n). The values would be recorded in the sequence

(x0,y0),   (x0+c,y0),   (x0+2c,y0), ...,   (x0+(n-1)*c,y0)
(x0,y0+c), (x0+c,y0+c), (x0+2c,y0+c), ..., (x0+(n-1)*c,y0+c)
...
(x0,y0+(m-1)*c), (x0+c, y0+(m-1)*c), ..., (x0+(n-1)*c,y0+(m-1)*c)

This sequence can be addressed with a single index k = 0, 1, ..., n*m-1.

Given an arbitrary point (x,y), you use simple calculations to find the indexes where its four nearest neighbors are stored: they lie between rows i and i+1 and columns j and j+1 where

i = Floor((y-y0)/c), j = Floor((x-x0)/c).

The corresponding offsets are ni + j, n(i+1) + j, ni + (j+1), and n(i+1) + (j+1).

This assumes x0 <= x < x+nc and y0 <= y < y+mc, which is worth checking in a general-purpose program. The coordinates of these four points are readily seen to be equal to

(x0 + j*c, y0 + i*c), (x0 + j*c, y0 + (i+1*c), etc.

Now that you can solve (a) efficiently, it remains to interpolate the values found, which I will denote z00, z01, z10, and z11, situated as suggested in this crude picture (and shown more precisely in the figure at the end):

z01    z11

     z
z00    z10

You wish to estimate the value of z. It is convenient to use bilinear interpolation. The formulas amount to the following: first interpolate linearly along the top and bottom edges of the square. Here, only the x coordinate is changing (from x0+j*c to x0+(j+1)*c). Rescale this to t varying from 0 to 1, where

t = (x - (x0+j*c))/c = (x-x0)/c - j.

The interpolated values are

z0 = (1-t)*z00 + t*z10, z1 = (1-t)*z01 + t*z11

Now interpolate vertically between z0 and z1. First rescale y to s varying from 0 to 1, where

s = (y-y0)/c - i

and find the interpolated value

z = (1-s)*z0 + s*z1.

Example

Let (x0,y0) = (-180,-90), c = 1/24 = 0.041667, n = 24*360 = 8640, and (x,y) = (54.4786674627, 17.0470721369) (that is, near latitude 17 degrees north and longitude 54.5 degrees east). Then

i = Floor((y-y0)/c) = Floor(24*(90 + 17.047...)) = 2569
j = Floor((x-x0)/c) = Floor(24*(180 + 54.478....)) = 5627.

The undulation for (i,j), namely z00, is found at offset 8640*i+j = 48,619,849. The undulation for (i+1,j), namely z01, is found at offset 8640*(i+1)+j = 48,628,489. The undulation for (i,j+1), namely z10, is found at offset 48,619,850, and the undulation for (i+1,j+1), namely z11, is found at offset 48,628,490. Suppose these have the values given; that is,

z00 = 31.945, z01 = 31.993, z10 = 31.866, and z11 = 31.911

To interpolate, compute

t = (x-x0)/c - j = 0.488019
s = (y-y0)/c - i = 0.129731

Therefore

z0 = 31.96842 and z1 = 31.88796

Finally

z = 31.958.