Ortho projection does not center correctly

obspy/core/tests/data/events_longitude_wrap.zmap

@@ -0,0 +1,15 @@
+-176.504900	-27.946500	2016.531543006224	7	13	5.700000	10.000000	13	4	25.56
+-176.550200	-28.145000	2016.531455983417	7	13	5.700000	10.000000	12	18	33.69
+-176.386100	-28.032700	2016.531442052153	7	13	6.300000	12.410000	12	11	13.15
+-176.319700	-28.221600	2016.531415032698	7	13	5.700000	17.180000	11	56	58.73
+-173.213900	-15.074500	2016.523422669057	7	10	5.800000	10.000000	13	44	41.01
+-172.855200	-15.023100	2016.523417376923	7	10	6.000000	10.000000	13	41	53.66
+-142.133600	-56.801600	2016.511046523667	7	6	5.900000	10.000000	1	1	57.59
+167.421100	-16.092800	2016.495845820368	6	30	5.800000	39.650000	11	30	34.87
+166.437400	-14.692200	2016.491664358809	6	28	5.800000	10.000000	22	46	47.02
+123.369100	23.689300	2016.477811002011	6	23	5.800000	10.000000	21	5	30.63
+151.853800	-3.462000	2016.471903727421	6	21	6.300000	365.310000	17	12	8.43
+168.759500	-20.207200	2016.467651261131	6	20	6.000000	15.000000	3	50	55.24
+169.075200	-20.280200	2016.465595388079	6	19	6.300000	13.000000	9	47	23.6
+166.341200	-11.599800	2016.461456986819	6	17	5.700000	43.250000	21	26	17.42
+168.829100	-18.762600	2016.452393323087	6	14	6.200000	111.000000	13	49	22.62

obspy/core/tests/test_event.py

@@ -501,6 +501,20 @@ def test_catalog_plot_ortho(self):
                      resolution='c', water_fill_color='#98b7e2', label=None,
                      color='date')
 
+    def test_catalog_plot_ortho_longitude_wrap(self):
+        """
+        Tests the catalog preview plot, ortho projection, some non-default
+        parameters, using Basemap, with longitudes that need the mean to be
+        computed in a circular fashion.
+        """
+        cat = read_events('/path/to/events_longitude_wrap.zmap', format='ZMAP')
+        with ImageComparison(self.image_dir,
+                             'catalog-basemap_long-wrap.png') as ic:
+            rcParams['savefig.dpi'] = 40
+            cat.plot(method='basemap', outfile=ic.name, projection='ortho',
+                     resolution='c', label=None, title='', colorbar=False,
+                     water_fill_color='b')
+
     def test_catalog_plot_local(self):
         """
         Tests the catalog preview plot, local projection, some more non-default
@@ -555,6 +569,20 @@ def test_catalog_plot_ortho(self):
                      resolution='c', water_fill_color='#98b7e2', label=None,
                      color='date')
 
+    def test_catalog_plot_ortho_longitude_wrap(self):
+        """
+        Tests the catalog preview plot, ortho projection, some non-default
+        parameters, using Cartopy, with longitudes that need the mean to be
+        computed in a circular fashion.
+        """
+        cat = read_events('/path/to/events_longitude_wrap.zmap', format='ZMAP')
+        with ImageComparison(self.image_dir,
+                             'catalog-cartopy_long-wrap.png') as ic:
+            rcParams['savefig.dpi'] = 40
+            cat.plot(method='cartopy', outfile=ic.name, projection='ortho',
+                     resolution='c', label=None, title='', colorbar=False,
+                     water_fill_color='b')
+
     def test_catalog_plot_local(self):
         """
         Tests the catalog preview plot, local projection, some more non-default

obspy/geodetics/base.py

@@ -16,6 +16,7 @@
 import warnings
 
 import numpy as np
+from scipy.stats import circmean
 
 from obspy.core.util.decorator import deprecated
 from obspy.core.util.misc import to_int_or_zero
@@ -358,6 +359,24 @@ def locations2degrees(lat1, long1, lat2, long2):
     return gd
 
 
+def mean_longitude(longitudes):
+    """
+    Compute sample mean longitude, assuming longitude in degrees from -180 to
+    180.
+
+    >>> lons = (-170.5, -178.3, 166)
+    >>> np.mean(lons)  # doctest: +SKIP
+    -60.933
+    >>> mean_longitude(lons)  # doctest: +ELLIPSIS
+    179.08509...
+
+    :type longitudes: :class:`~numpy.ndarray` (or list, ..)
+    :param longitudes: Geographical longitude values ranging from -180 to 180
+        in degrees.
+    """
+    return circmean(np.array(longitudes), low=-180, high=180)
+
+
 if __name__ == '__main__':
     import doctest
     doctest.testmod(exclude_empty=True)

obspy/imaging/maps.py

@@ -27,6 +27,7 @@
 
 from obspy import UTCDateTime
 from obspy.core.util.base import get_basemap_version, get_cartopy_version
+from obspy.geodetics.base import mean_longitude
 
 BASEMAP_VERSION = get_basemap_version()
 if BASEMAP_VERSION:
@@ -316,7 +317,7 @@ def _plot_basemap_into_axes(
             bmap = Basemap(projection='ortho',
                            resolution=_BASEMAP_RESOLUTIONS[resolution],
                            area_thresh=1000.0, lat_0=round(np.mean(lats), 4),
-                           lon_0=round(np.mean(lons), 4), ax=ax)
+                           lon_0=round(mean_longitude(lons), 4), ax=ax)
         elif projection == 'local':
             if min(lons) < -150 and max(lons) > 150:
                 max_lons = max(np.array(lons) % 360)
@@ -576,7 +577,7 @@ def plot_cartopy(lons, lats, size, color, labels=None, projection='global',
         proj = ccrs.Mollweide(**proj_kwargs)
     elif projection == 'ortho':
         proj_kwargs['central_latitude'] = np.mean(lats)
-        proj_kwargs['central_longitude'] = np.mean(lons)
+        proj_kwargs['central_longitude'] = mean_longitude(lons)
         proj = ccrs.Orthographic(**proj_kwargs)
     elif projection == 'local':
         if min(lons) < -150 and max(lons) > 150:
@@ -619,7 +620,7 @@ def plot_cartopy(lons, lats, size, color, labels=None, projection='global',
     elif isinstance(projection, type):
         if 'central_longitude' in proj_kwargs:
             if proj_kwargs['central_longitude'] == 'auto':
-                proj_kwargs['central_longitude'] = np.mean(lons)
+                proj_kwargs['central_longitude'] = mean_longitude(lons)
         if 'central_latitude' in proj_kwargs:
             if proj_kwargs['central_latitude'] == 'auto':
                 proj_kwargs['central_latitude'] = np.mean(lats)