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:
Best Answer
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
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
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
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):
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
The interpolated values are
Now interpolate vertically between z0 and z1. First rescale y to s varying from 0 to 1, where
and find the interpolated value
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
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,
To interpolate, compute
Therefore
Finally