Skip to content

Commit

Permalink
Replace geographiclib with pyproj for performing lat/lon calculations…
Browse files Browse the repository at this point in the history
… -- pyproj is much much faster, and already brought in by Basemap
  • Loading branch information
adam-iris committed Oct 19, 2017
1 parent 59d34b7 commit b7c03ad
Show file tree
Hide file tree
Showing 4 changed files with 27 additions and 15 deletions.
3 changes: 1 addition & 2 deletions environment.yml
Original file line number Diff line number Diff line change
Expand Up @@ -14,5 +14,4 @@ dependencies:
- pyqt=4.11.4
- qtconsole
- pyzmq
- pip:
- geographiclib
- pyproj
5 changes: 4 additions & 1 deletion installer/anaconda/meta.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,11 @@ requirements:
run:
- python
- obspy
- pyqt 4.11.4
- pyqt 4.11*
- qtconsole
- basemap
- pyproj
- pillow
- mock # [py2k]

test:
Expand Down
2 changes: 1 addition & 1 deletion installer/environment.yml
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,6 @@ dependencies:
- pyqt=4.11.4
- qtconsole
- pyzmq
- pyproj
- pip:
- pyweed==0.4.dev4
- geographiclib
32 changes: 21 additions & 11 deletions pyweed/pyweed_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,14 +15,18 @@
import os
import logging
import re
from geographiclib.geodesic import Geodesic
from pyproj import Geod
from obspy.taup.tau import TauPyModel
from future.moves.urllib.parse import urlencode

LOGGER = logging.getLogger(__name__)
GEOD = Geodesic.WGS84
GEOD = Geod(ellps='WGS84')
TAUP = TauPyModel()

# Rough meters/degree calculation
M_PER_DEG = (GEOD.inv(0, 0, 0, 1)[2] + GEOD.inv(0, 0, 1, 0)[2]) / 2



class OutputFormat(object):
"""
Expand Down Expand Up @@ -278,8 +282,9 @@ def get_distance(lat1, lon1, lat2, lon2):
"""
Get the distance between two points in degrees
"""
dist = GEOD.Inverse(lat1, lon1, lat2, lon2, Geodesic.DISTANCE)
return dist['a12']
# NOTE that GEOD takes longitude first!
dist = GEOD.inv(lon1, lat1, lon2, lat2)
return dist[2]


def get_arrivals(distance, event_depth):
Expand Down Expand Up @@ -311,13 +316,18 @@ def get_bounding_circle(lat, lon, radius, num_points=36):
"""
Returns groups of lat/lon pairs representing a circle on the map
"""
return [
(g['lat2'], g['lon2'])
for g in [
GEOD.ArcDirect(lat, lon, (i * 360) / num_points, radius)
for i in range(0, num_points + 1)
]
]
radius_meters = radius * M_PER_DEG
# NOTE that GEOD takes longitude first!
trans = GEOD.fwd(
[lon] * num_points,
[lat] * num_points,
list(((i * 360) / num_points) for i in range(num_points)),
[radius_meters] * num_points
)
points = list(zip(trans[1], trans[0]))
# We need to complete the circle by adding the first point again as the last point
points.append(points[0])
return points


def get_service_url(client, service, parameters):
Expand Down

0 comments on commit b7c03ad

Please sign in to comment.