USGS Earthquake GPS Metadata: Buffer Zones & Python

Question

I'm making a python web crawler that gets USGS earthquake GPS data, and I'm wondering how I should process the latitude longitude pairs to see if it falls within a specified radius of a specified GPS location. What is the mathematical model---pseudocode really---that needs to be used for this problem?

Answer 1

If you don't want to write your own code you can use the GDAL Python bindings to do that.

Let's say this is your GPS point with lat and lon:

from osgeo import ogr, osr

# your gps earthquacke point
gps_lat = 47.86
gps_lon = 12.66 
gps_point = ogr.Geometry(ogr.wkbPoint)
gps_point.AddPoint(gps_lon,gps_lat)

So your earthquake happens somewhere in the Alps

Now we create a second point, the center of your buffer

# your location that is the center of the buffer
location_lat = 48.140432
location_lon = 11.615295 
location_point = ogr.Geometry(ogr.wkbPoint)
location_point.AddPoint(location_lon,location_lat)

The center of the buffer will be Munich.

In order to create a buffer in meters around Munich we need to transform our points to a projected coordinated system in meters

# create a transformation from EPSG:4326 (lat/long) to EPSG:3035 (Projected coordinate system for Europe in meters)
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(4326)
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(3035)
coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# transform points to to EPSG:3035
location_point.Transform(coordTransform)
gps_point.Transform(coordTransform)

Now we can create a buffer around Munich in meters. Lets say 200km:

# buffer location to create buffer circle
buffer_distance = 200*1000 # in meter
buffer_polygon = location_point.Buffer(buffer_distance)

Last we can check if the earthquake point is in the buffer:

print buffer_polygon.Contains(gps_point)
>>> True