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?
-
You've drawn it as a plane. Did you want a plane or great circle. Assuming that latter: Haversine (en.wikipedia.org/wiki/Haversine_formula) – BradHards Apr 16 '14 at 6:27
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