Python, numpy, scipy: How to exclude a location error from GPS coordinates? (averaging lon, lat)

Question

I use Python:

I have 2 arrays of GPS points - lon and lat (more than 500 000 points).

I have 1 array of date-time.

lon = numpy.array(lon)
lat = numpy.array(lat)
dt = numpy.array(dt)

I have a location error (GPS sensor error). For example 15 meters.

GPS_sensor_error = 0.015    

I need to exclude a GPS_sensor_error from coordinates that there were no asterisks on a track.

(I don't draw a point with identical coordinates)

How I am able to do it?

Now:

1. I calculate distance between points.

2. I find the minimum distance, if it less GPS_sensor_error, then I average lon, lat.

3. repeat 1.

4. repeat 2.

5. repeat until all distances won't be more GPS_sensor_error

Update:

     lon = numpy.array()
     lat = numpy.array()

     flag = True
     while flag:
        lon1 = lon[:-1]
        lon2 = lon[1:]
        lat1 = lat[:-1]
        lat2 = lat[1:]

        '''distance'''
        x = (lon2 - lon1)
        y = (lat2 - lat1)
        d = numpy.sqrt(x * x + y * y)

        min = numpy.min(d)
        if min < GPS_sensor_error:
            j = numpy.where(d == min)[0][0]

            lon[j] = (lon[j] + lon[j + 1]) / 2
            lat[j] = (lat[j] + lat[j + 1]) / 2

            lon = numpy.delete(lon, j + 1)
            lat = numpy.delete(lat, j + 1)

        else:
            flag = False

Bypass on all points works at a pure python very long... Prompt please, how to implement it using scipy, numpy?

Thanks

P.s. probably already there is a GPS filter in scipy, numpy?

Answer 1

From a data science perspective what you are doing is not correct. You cant just use the average error distance as a cutoff and think your data will be more correct. The two points you are comparing can have an error more or less than 15 m they can shift toward each other or move away of each other. And if you don't have another exact dataset there is no way of telling what would be the correct point. You can't make this dataset more precise.

However I think you goal is to simplify your dataset, not to make it more accurate. For that you can use the Douglas–Peucker algorithm. I would suggest you load your data in an Postgis enabled database (Postgresql + postgis) and then use the simplify function. This will require some db setup time, but then it will speed you up greatly. However if you want it in pure python this SO question has a very nice snippet.

BTW. If your are doing distance calculations with lat,lon do not use Pythagoras. It is not valid since lat,lon is not Euclidean. Use the haversine algorithm.

Answer 2

You can easily do all your calculations using numpy primitives only and no python looping.

First define your distance function as a function that operates on numpy arrays (I assume you did that already ..):

def dist(lon1, lat1, lon2, lat2):
    """Compute the distance between (lon1, lat1) and (lon2, lat2). 
       Both may be numpy arrays."""
    ...

Then apply it to your data as so:

d = dist(lon[:-1], lat[:-1], lon[1:], lat[1:])

This notation means you will compare the i^th point to the i+1^th.

Next find the indices where d is greater than your threshold:

I = d > GPS_sensor_error

Now keep only those and the first point!

lon_out = numpy.hstack([[lon[0]], lon[1:][I]]) # could also use numpy.where
lat_out = numpy.hstack([[lat[0]], lat[1:][I]])

Update:

If you want to keep the same number of points, ie set lon[i] to the last good value, use the following trick instead of the previous two lines:

idx, = numpy.where(I)
idx = numpy.hstack([[0], idx])
J = numpy.cumsum(I) # the trick
lon_out = lon[idx[J]]
lat_out = lat[idx[J]]