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:
if not x1 < x < x2 or not y1 < y < y2
to<=
, you should then also confirm that x2 != x1 and y2 != y1.