Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
1  
Thank you for posting the code. Please note that it will fail to interpolate whenever the point (x,y) falls on the boundary of a grid cell: that's an important flaw that should be fixed; such circumstances do occur in practice. If you try the quick fix of changing the strict inequalities in the test if not x1 < x < x2 or not y1 < y < y2 to <=, you should then also confirm that x2 != x1 and y2 != y1. –  whuber Dec 29 '11 at 5:30
    
@whuber Thanks for vigilant comment. –  daikini Dec 29 '11 at 13:48

1 Answer 1

up vote 8 down vote accepted

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 n*i + j, n*(i+1) + j, n*i + (j+1), and n*(i+1) + (j+1).

This assumes x0 <= x < x+n*c and y0 <= y < y+m*c, 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.

Figure

share|improve this answer
2  
thank you very much for your complete lecture and your time. –  daikini Oct 17 '11 at 15:25
    
how the coordinates of origin of my geoid grid determine accuracy of computation? –  daikini Oct 17 '11 at 19:51
    
@user Inaccuracy in the origin translates directly to the same inaccuracy in the coordinates (x,y) of the GPS point. How much of a difference that makes depends on how variable the surrounding values z00, ..., z11 are. But usually there is no effective inaccuracy in grid coordinates: the grid is established, measurements are taken at the grid nodes, and any inaccuracy due to imprecise spatial location in the field is just rolled up into the overall inaccuracy of the measurements themselves. –  whuber Oct 17 '11 at 20:32
2  
For the record, check also geographiclib.sourceforge.net/cgi-bin/GeoidEval –  markusN Dec 30 '11 at 14:45

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.