skymap.py 24.7 KB
Newer Older
DavidWalz's avatar
DavidWalz committed
1
"""
2
Module to visualize arrival directions or heatmaps on the sphere.
DavidWalz's avatar
DavidWalz committed
3
"""
4

5
import healpy as hp
6
7
import matplotlib.pyplot as plt
import numpy as np
DavidWalz's avatar
skymap  
DavidWalz committed
8

9
from astrotools import coord, healpytools as hpt
DavidWalz's avatar
skymap  
DavidWalz committed
10

MarCus's avatar
MarCus committed
11

12
def scatter(v, c=None, cblabel='log$_{10}$(Energy / eV)', opath=None, fig=None, **kwargs):
DavidWalz's avatar
skymap  
DavidWalz committed
13
    """
14
    Scatter plot of events with arrival directions x,y,z and colorcoded energies.
15

16
    :param v: array of shape (3, n) pointing into directions of the events
17
18
19
    :param c: quantity that is supposed to occur in colorbar, e.g. energy of the cosmic rays
    :param cblabel: colorbar label
    :param opath: if not None, saves the figure to the given opath (no returns)
20
    :param fig: figure to plot in, creates new figure if None
Marcus Wirtz's avatar
Marcus Wirtz committed
21
    :param kwargs: additional named keyword arguments
22

23
           - figsize: figure size as input for plt.figure()
24
           - cmap: colormap
25
           - cbar: if True includes a colobar
26
           - cticks: sets ticks of colormap
27
28
           - mask_alpha: alpha value for maskcolor
           - fontsize: scale the general fontsize
Marcus Wirtz's avatar
Marcus Wirtz committed
29
30
31
32
33
           - dark_grid: if True paints a dark grid (useful for bright maps)
           - gridcolor: Color of the grid.
           - gridalpha: Transparency value of the gridcolor.
           - tickcolor: Color of the ticks.
           - tickalpha: Transparency of the longitude ticks.
34
35
36
           - plane: plots 'SGP' or 'GP' or both (list) into plot
           - planecolor: color of plane
           - coord_system: default galactic ('gal') / equatorial ('eq')
Marcus Wirtz's avatar
Marcus Wirtz committed
37
    :return: figure, axis of the scatter plot
DavidWalz's avatar
skymap  
DavidWalz committed
38
    """
39

DavidWalz's avatar
DavidWalz committed
40
    lons, lats = coord.vec2ang(v)
41

42
    fontsize = kwargs.pop('fontsize', 26)
43
    kwargs.setdefault('s', 8)
44
45
    if 'marker' not in kwargs:
        kwargs.setdefault('lw', 0)
46
    cbar = kwargs.pop('cbar', True) and isinstance(c, (list, tuple, np.ndarray))
47
    if cbar:
48
49
        vmin = kwargs.pop('vmin', smart_round(np.min(c[np.isfinite(c)]), upper_border=False))
        vmax = kwargs.pop('vmax', smart_round(np.max(c[np.isfinite(c)]), upper_border=True))
50
        step = smart_round((vmax - vmin) / 5., order=1)
51
        cticks = kwargs.pop('cticks', np.round(np.arange(vmin, vmax, step), int(np.round(-np.log10(step), 0))))
52
        clabels = kwargs.pop('clabels', cticks)
53

Marcus Wirtz's avatar
Marcus Wirtz committed
54
55
56
57
58
59
    # read keyword arguments for the grid
    dark_grid = kwargs.pop('dark_grid', True)
    gridcolor = kwargs.pop('gridcolor', 'lightgray' if dark_grid is None else 'black')
    gridalpha = kwargs.pop('gridalpha', 0.5 if dark_grid is None else 0.4)
    tickcolor = kwargs.pop('tickcolor', 'lightgray' if dark_grid is None else 'black')
    tickalpha = kwargs.pop('tickalpha', 0.5 if dark_grid is None else 1)
60
    planecolor = kwargs.pop('planecolor', 'darkgray')
61
62
63
64
    plane = kwargs.pop('plane', None)
    coord_system = kwargs.pop('coord_system', 'gal')

    if coord_system == 'eq':
65
        lons, lats = coord.vec2ang(coord.gal2eq(coord.ang2vec(lons, lats)))
Teresa Bister's avatar
Teresa Bister committed
66

67
68
    # mimic astronomy convention: positive longitudes evolving to the left with respect to GC
    lons = -lons
Marcus Wirtz's avatar
Marcus Wirtz committed
69

70
    # plot the events
71
    fig = plt.figure(figsize=kwargs.pop('figsize', [12, 6])) if fig is None else fig
72
73
74
    ax = fig.add_axes([0.1, 0.1, 0.85, 0.9], projection="hammer")
    events = ax.scatter(lons, lats, c=c, **kwargs)

75
    if cbar:
76
77
        cbar = plt.colorbar(events, orientation='horizontal', shrink=0.85, pad=0.05,
                            aspect=30, cmap=kwargs.get('cmap'), ticks=cticks)
78
        cbar.set_label(cblabel, fontsize=fontsize)
79
        events.set_clim(vmin, vmax)
80
        cbar.ax.tick_params(labelsize=fontsize - 4)
81
        cbar.set_ticklabels(clabels)
82
        cbar.draw_all()
83

Marcus Wirtz's avatar
Marcus Wirtz committed
84
85
86
87
    # Setup the grid
    plt.xticks(fontsize=fontsize)
    plt.yticks(fontsize=fontsize)
    plot_grid(gridcolor=gridcolor, gridalpha=gridalpha, tickalpha=tickalpha, tickcolor=tickcolor, fontsize=fontsize)
88

89
90
    if plane is not None:
        plot_plane(planecolor, coord_system, plane)
91

marcus's avatar
marcus committed
92
93
94
    if opath is not None:
        plt.savefig(opath, bbox_inches='tight')
        plt.clf()
DavidWalz's avatar
skymap  
DavidWalz committed
95

Marcus Wirtz's avatar
Marcus Wirtz committed
96
97
    return fig, ax

DavidWalz's avatar
skymap  
DavidWalz committed
98

99
100
101
102
103
104
105
def heatmap(m, opath=None, label='entries', mask=None, maskcolor='white', **kwargs):
    """
    Heatmap plot of binned data m. For exmaple usage see: cosmic_rays.plot_healpy_map()

    :param m: Array with size npix for an arbitrary healpy nside.
    :param opath: if not None, saves the figure to the given opath (no returns)
    :param label: label for the colormap
106
    :param mask: either boolean mask that paints certain pixels different or condition for m
107
    :param maskcolor: which color to paint the mask
108
    :param kwargs: additional named keyword arguments (non-listed keywords are passed to matplotlib.pcolormesh)
109

110
           - figsize: figure size as input for plt.figure()
111
           - cmap: colormap
112
113
114
           - vmin: Lower cut of the colorbar scale
           - vmax: Upper cut of the colorbar scale
           - cbticks: Ticks of the colorbar
115
116
           - mask_alpha: alpha value for maskcolor
           - fontsize: scale the general fontsize
117
           - xresolution: Scales the resolution of the plot (default: 800)
118
           - dark_grid: if True paints a dark grid (useful for bright maps)
Marcus Wirtz's avatar
Marcus Wirtz committed
119
120
121
122
           - gridcolor: Color of the grid.
           - gridalpha: Transparency value of the gridcolor.
           - tickcolor: Color of the ticks.
           - tickalpha: Transparency of the longitude ticks.
123
124
125
           - plane: plots 'SGP' or 'GP' or both (list) into plot
           - planecolor: color of plane
           - coord_system: default galactic ('gal') / equatorial ('eq')
126
127
128
    :return: figure of the heatmap, colorbar
    """

129
    # read general keyword arguments
Marcus Wirtz's avatar
Marcus Wirtz committed
130
    fontsize = kwargs.pop('fontsize', 26)
131

132
133
134
135
136
    # create the meshgrid and project the map to a rectangular matrix (xresolution x yresolution)
    xresolution = kwargs.pop('xresolution', 800)
    yresolution = xresolution // 2
    theta = np.linspace(np.pi, 0, yresolution)
    phi = np.linspace(-np.pi, np.pi, xresolution)
137
138
    longitude = np.deg2rad(np.linspace(-180, 180, xresolution))
    latitude = np.deg2rad(np.linspace(-90, 90, yresolution))
139
140
141
142
143

    phi_grid, theta_grid = np.meshgrid(phi, theta)
    grid_pix = hp.ang2pix(hp.get_nside(m), theta_grid, phi_grid)

    # management of the colormap
144
    cmap = kwargs.pop('cmap', 'viridis')
145
146
147
148
    coord_system = kwargs.pop('coord_system', 'gal')
    planecolor = kwargs.pop('planecolor', 'darkred')
    plane = kwargs.pop('plane', None)

149
150
151
152
153
154
    if isinstance(cmap, str):
        cmap = plt.cm.get_cmap(cmap)
    if mask is not None:
        if not hasattr(mask, 'size'):
            mask = m == mask
        m = np.ma.masked_where(mask, m)
155
        cmap.set_bad(maskcolor, alpha=kwargs.pop('mask_alpha', 1))
156
157
158
159
160
161
162

    if coord_system == 'eq':
        hrot = hp.Rotator(coord=['G', 'C'], inv=True)
        theta, phi = hp.pix2ang(hpt.get_nside(m), np.arange(len(m)))
        g0, g1 = hrot(theta, phi)
        pix0 = hp.ang2pix(hpt.get_nside(m), g0, g1)
        m = m[pix0]
163
    grid_map = m[grid_pix]
164

165
    finite = np.isfinite(m)
166
167
    vmin = kwargs.pop('vmin', smart_round(np.min(m[finite]), upper_border=False))
    vmax = kwargs.pop('vmax', smart_round(np.max(m[finite]), upper_border=True))
168
    cbticks = kwargs.pop('cbticks', [vmin, (vmin + vmax) / 2, vmax])
169
170

    # read keyword arguments for the grid
171
    dark_grid = kwargs.pop('dark_grid', None)
172
173
174
    gridcolor = kwargs.pop('gridcolor', 'lightgray' if dark_grid is None else 'black')
    gridalpha = kwargs.pop('gridalpha', 0.5 if dark_grid is None else 0.4)
    tickcolor = kwargs.pop('tickcolor', 'lightgray' if dark_grid is None else 'black')
175
    tickalpha = kwargs.pop('tickalpha', 0.5 if dark_grid is None else 1)
DavidWalz's avatar
DavidWalz committed
176

177
    # Start plotting
178
    fig = plt.figure(figsize=kwargs.pop('figsize', (12, 12)))
179
    fig.add_subplot(111, projection='hammer')
DavidWalz's avatar
DavidWalz committed
180
    # flip longitude to the astro convention
181
    # rasterized makes the map bitmap while the labels remain vectorial
182
    image = plt.pcolormesh(longitude[::-1], latitude, grid_map, rasterized=True, vmin=vmin,
Marcus Wirtz's avatar
Marcus Wirtz committed
183
                           vmax=vmax, cmap=cmap, edgecolor='face', **kwargs)
184
    cb = fig.colorbar(image, ticks=cbticks, orientation='horizontal', aspect=30, shrink=0.9, pad=0.05)
MarCus's avatar
MarCus committed
185
    cb.solids.set_edgecolor("face")
186
187
    cb.set_label(label, fontsize=fontsize)
    cb.ax.tick_params(axis='x', direction='in', size=3, labelsize=fontsize - 4)
188

189
190
    if plane is not None:
        plot_plane(planecolor, coord_system, plane)
191

MarCus's avatar
MarCus committed
192
    # Setup the grid
Marcus Wirtz's avatar
Marcus Wirtz committed
193
    plot_grid(gridcolor=gridcolor, gridalpha=gridalpha, tickalpha=tickalpha, tickcolor=tickcolor, fontsize=fontsize)
marcus's avatar
marcus committed
194
195
196
197

    if opath is not None:
        plt.savefig(opath, bbox_inches='tight')
        plt.clf()
198
199

    return fig, cb
200
201


Marcus Wirtz's avatar
Marcus Wirtz committed
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
def plot_grid(lon_ticks=None, lat_ticks=None, lon_grid=30, lat_grid=30, fontsize=26, **kwargs):
    """
    Plot a grid on the skymap.

    :param lon_ticks: Set the label ticks for the longitudes (default: [90, 0, -90]).
    :param lat_ticks: Set the label ticks for the latitude (default: [-60, -30, 0, 30, 60]).
    :param lon_grid: Distances between the grid lines in longitude in degrees (default: 30 deg).
    :param lat_grid: Distances between the grid lines in latitude in degrees (default: 30 deg).
    :param kwargs: additional named keyword arguments

            - gridcolor: Color of the grid.
            - gridalpha: Transparency value of the gridcolor.
            - tickcolor: Color of the ticks.
            - tickalpha: Transparency of the longitude ticks.
    """
    lon_ticks = [90, 0, -90] if lon_ticks is None else lon_ticks
    lat_ticks = [-60, -30, 0, 30, 60] if lat_ticks is None else lat_ticks
    plt.gca().set_longitude_grid(lon_grid)
    plt.gca().set_latitude_grid(lat_grid)
    plt.gca().set_longitude_grid_ends(89)

    plt.grid(alpha=kwargs.pop('gridalpha', 0.5), color=kwargs.pop('gridcolor', 'lightgray'))
    plt.gca().set_xticklabels(np.hstack([(r'', r'', r'%d$^{\circ}$' % lon) for lon in lon_ticks]),
                              alpha=kwargs.pop('tickalpha', 0.5), fontsize=fontsize)
    plt.gca().tick_params(axis='x', colors=kwargs.pop('tickcolor', 'lightgray'))
    plt.gca().set_yticklabels([r'%d$^{\circ}$' % lat for lat in lat_ticks], fontsize=fontsize)


230
def plot_plane(planecolor=0.5, coord_system='gal', plane='SGP', **kwargs):
231
    """
232
    Plot a the supergalactic plane onto skymap.
233
234

    :param planecolor: color of plane
235
236
237
    :param coord_system: default galactic ('gal') / equatorial ('eq')
    :param plane: plots 'SGP' or 'GAL' or both (list) into plot
    :param kwargs: additional named keyword arguments passed to plt.plot()
238
239
240
241
242
    """
    phi0 = np.linspace(0, 2*np.pi, 100)
    if coord_system.upper() == 'GAL':
        # only plotting the SGP makes sense
        phi, theta = coord.vec2ang(coord.sgal2gal(coord.ang2vec(phi0, np.zeros_like(phi0))))
243
244
        kwargs.setdefault('color', planecolor)
        plt.plot(-np.sort(phi), theta[np.argsort(phi)], **kwargs)
245
246
247
248

    elif coord_system.upper() == 'EQ':
        if 'SGP' in plane:
            phi, theta = coord.vec2ang(coord.gal2eq(coord.sgal2gal(coord.ang2vec(phi0, np.zeros_like(phi0)))))
249
250
            kwargs.setdefault('color', planecolor)
            plt.plot(-np.sort(phi), theta[np.argsort(phi)], **kwargs)
251
252
        if 'GAL' in plane:
            phi, theta = coord.vec2ang(coord.gal2eq(coord.ang2vec(phi0, np.zeros_like(phi0))))
253
254
            kwargs.setdefault('color', '0.5')
            plt.plot(-np.sort(phi), theta[np.argsort(phi)], **kwargs)
255
256
257
258
259
260
        else:
            raise Exception("plane type not understood, use GP or SGP or list!")
    else:
        raise Exception("coord system not understood, use eq or gal!")


261
def plot_tissot(vec_c, r, res=1e-2, **kwargs):
262
263
    """
    Plot a circle onto skymap.
264

265
266
267
268
269
    :param vec_c: vector pointing to the center of the circle
    :param r: radius of the circle
    :param res: resolution of the circle approximation (in radian)
    :param kwargs: additional named keyword arguments passed to plt.plot()
    """
270
271
272
273
274
275
276
    lon, lat = coord.vec2ang(vec_c)
    vec_ref = coord.rotate(vec_c, coord.sph_unit_vectors(lon, lat)[2], r)
    psis = np.arange(0, 2*np.pi, res)
    lons, lats = coord.vec2ang(coord.rotate(vec_ref, vec_c, psis))
    plt.plot(-lons, lats, **kwargs)


277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
def smart_round(v, order=2, upper_border=True):
    """
    Rounds a value v such that it can be used e.g. for colorbars

    :param v: scalar value which should be rounded
    :type v: int, float
    :param upper_border: round such that the value can be used as an upper border of an interval, default=True
    :param order: number of digits to round to, default=2
    :return: rounded value
    :rtype: int, float

    This function has been tested on the following numbers (with all upper_border presented here):

    .. code-block:: python

        :linenos:
293
        >> from astrotools.skymap import smart_round
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
        >> smart_round(100000), smart_round(100000, upper_border=False)
        100000.0, 100000.0
        >> smart_round(100001), smart_round(100001, upper_border=False)
        101000.0, 100000.0
        >> smart_round(-100001), smart_round(-100001, upper_border=False)
        -100000.0, -100000.0
        >> smart_round(2.23), smart_round(2.23, upper_border=False)
        2.23, 2.23
        >> smart_round(2.230), smart_round(2.230, upper_border=False)
        2.23, 2.23
        >> smart_round(2.231), smart_round(2.231, upper_border=False)
        2.24, 2.23
        >> smart_round(-2.230), smart_round(-2.230, upper_border=False)
        -2.23, -2.23
        >> smart_round(-2.231), smart_round(-2.231, upper_border=False)
        -2.23, -2.24
        >> smart_round(0.930001), smart_round(0.930001, upper_border=False)
        0.94, 0.93
        >> smart_round(-0.930001), smart_round(-0.930001, upper_border=False)
        -0.93, -0.94
    """
    if v == 0:
        return 0
    o = np.log10(np.fabs(v))
    f = 10 ** (-int(o) + order)
    if upper_border:
        return np.ceil(v * f) / f
    return np.floor(v * f) / f


324
def skymap(m, **kwargs):
Lukas Geiger's avatar
Lukas Geiger committed
325
    """ Deprecated funcion -> See heatmap() """
Marcus Wirtz's avatar
Marcus Wirtz committed
326
327
    print("User warning: function skymap() is deprecated and will stop existing in v2.0.0. \
            Please use heatmap() in future!")
328
    return heatmap(m, **kwargs)
Marcus Wirtz's avatar
Marcus Wirtz committed
329
330
331


def healpy_map(m, **kwargs):
332
    """ Forwards to function heatmap() """
Marcus Wirtz's avatar
Marcus Wirtz committed
333
334
335
336
    return heatmap(m, **kwargs)


def eventmap(v, **kwargs):
337
    """ Forwards to function scatter() """
Marcus Wirtz's avatar
Marcus Wirtz committed
338
    return scatter(v, **kwargs)
339
340
341
342


class PlotSkyPatch:

343
344
345
346
    """
    Class to plot a close-up look of a region of interest (ROI) in the sky.
    To use this class you need to install the Basemap package:
    https://matplotlib.org/basemap/users/installing.html
347

348
    .. code-block:: python
349

350
351
352
353
354
355
356
357
        from astrotools.skymap import PlotSkyPatch
        patch = PlotSkyPatch(lon0, lat0, r_roi, title='My Skypatch')
        mappable = patch.plot_crs("/path/to/cosmic_rays.CosmicRaysSets.npz", set_idx=0)
        patch.mark_roi()
        patch.plot_grid()
        patch.colorbar(mappable)
        patch.savefig("/tmp/test-skypatch.png")
    """
358

359
360
    def __init__(self, lon_roi, lat_roi, r_roi, ax=None, title=None, **kwargs):
        """
361
        :param lon_roi: Longitude of center of ROI in radians (0..2*pi)
362
        :param lat_roi: Latitude of center of ROI in radians (0..2*pi)
363
364
365
366
367
        :param r_roi: Radius of ROI to be plotted (in radians)
        :param ax: Matplotlib axes in case you want to plot on certain axes
        :param title: Optional title of plot (plotted in upper left corner)
        :param kwargs: keywords passed to matplotlib.figure()
        """
368
        from mpl_toolkits.basemap import Basemap    # pylint: disable=import-error,no-name-in-module
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
        import matplotlib as mpl

        with_latex_style = {
            "text.usetex": True,
            "font.family": "serif",
            "axes.labelsize": 30,
            "font.size": 30,
            "legend.fontsize": 30,
            "xtick.labelsize": 26,
            "ytick.labelsize": 26,
            "legend.fancybox": False,
            "lines.linewidth": 3.0,
            "patch.linewidth": 3.0}

        mpl.rcParams.update(with_latex_style)

385
386
387
        assert (isinstance(lon_roi, (float, int))) and (isinstance(lat_roi, (float, int))) and \
            (isinstance(r_roi, (float, int))), "Keywords 'lon_roi', 'lat_roi' and 'r_roi' have to be floats or ints!"

388
389
390
391
392
393
394
395
        self.vec_0 = coord.ang2vec(lon_roi, lat_roi)
        self.lon_0 = np.rad2deg(lon_roi)
        self.lat_0 = np.rad2deg(lat_roi)
        self.r_roi = r_roi

        self.scale = 5500000 * (r_roi / 0.3)

        self.fig = None
396
        self.ax = ax
397
398
399
400
401
402
403
        if ax is None:
            kwargs.setdefault('figsize', [8, 8])
            self.fig = plt.figure(**kwargs)
            self.ax = plt.axes()

        self.title = title
        if title is not None:
404
            self.text(0.02, 0.98, title, verticalalignment='top', fontsize=36)
405
406
407
408
409
410
411
412

        self.m = Basemap(width=self.scale, height=self.scale, resolution='l', projection='stere', celestial=True,
                         lat_0=self.lat_0, lon_0=-360 - self.lon_0 if self.lon_0 < 0 else -self.lon_0, ax=ax)

    def plot_crs(self, crs, set_idx=0, zorder=0, cmap='viridis', **kwargs):
        """
        Plot cosmic ray events in the sky.

413
        :param crs: Either cosmic_rays.CosmicRaysBase or cosmic_rays.CosmicRaysSets object (or path) or dict object
414
415
416
417
418
        :param set_idx: In case of CosmicRaysSets object, chose the respective set index
        :param zorder: Usual matplotlib zorder keyword (order of plotting)
        :param cmap: Matplotlib colormap object or string
        """
        if isinstance(crs, str):
419
            from astrotools import cosmic_rays
420
421
422
423
            try:
                crs = cosmic_rays.CosmicRaysBase(crs)
            except AttributeError:
                crs = cosmic_rays.CosmicRaysSets(crs)
424
425

        if hasattr(crs, 'type') and (crs.type == "CosmicRaysSet"):
426
427
            crs = crs[set_idx]

428
        if 'log10e' in crs.keys():
429
            log10e = crs['log10e']
430
            assert np.all(log10e < 25), "Input energies ('log10e' key) are too high for being plotted"
431
432
            kwargs.setdefault('s', 10**(log10e - 18.))
            kwargs.setdefault('c', log10e)
433
434
        kwargs.setdefault('lw', 0)

435
        return self.scatter(crs['lon'], crs['lat'], zorder=zorder, cmap=cmap, **kwargs)
436
437
438
439
440
441
442
443
444
445
446
447
448

    def plot(self, lons, lats, **kwargs):
        """ Replaces matplotlib.pyplot.plot() function """
        kwargs.setdefault('rasterized', True)
        x, y = self.m(np.rad2deg(lons), np.rad2deg(lats))
        return self.m.plot(x, y, **kwargs)

    def scatter(self, lons, lats, **kwargs):
        """ Replaces matplotlib.pyplot.scatter() function """
        kwargs.setdefault('rasterized', True)
        x, y = self.m(np.rad2deg(lons), np.rad2deg(lats))
        return self.m.scatter(x, y, **kwargs)

449
    def tissot(self, lon, lat, radius, npts=1000, **kwargs):
450
451
452
453
        """ Replaces the Basemap tissot() function (plot circles) """
        kwargs.setdefault('fill', False)
        kwargs.setdefault('lw', 1)
        kwargs.setdefault('color', 'grey')
454
        return self.m.tissot(np.rad2deg(lon), np.rad2deg(lat), np.rad2deg(radius), npts, **kwargs)
455

456
    def mark_roi(self, alpha=0.4, **kwargs):
457
458
459
460
461
462
463
464
465
        """
        Marks the ROI by a circle ans shades cosmic rays outside the ROI.

        :param kwargs: Passed to Basemaps tissot() function
        """
        from matplotlib import path, collections
        kwargs.setdefault('lw', 2)
        kwargs.setdefault('zorder', 3)
        try:
466
            t = self.tissot(np.deg2rad(self.lon_0), np.deg2rad(self.lat_0), self.r_roi, **kwargs)
467
468
469
470
471
472
473
474
            xyb = np.array([[0., 0.], [1., 0.], [1., 1.], [0., 1.], [0., 0.]]) * self.scale

            p = path.Path(np.concatenate([xyb, t.get_xy()[::-1]]))
            p.codes = np.ones(len(p.vertices), dtype=p.code_type) * p.LINETO
            p.codes[0] = path.Path.MOVETO
            p.codes[4] = path.Path.CLOSEPOLY
            p.codes[5] = path.Path.MOVETO
            p.codes[-1] = path.Path.CLOSEPOLY
475
            col = collections.PathCollection([p], facecolor='white', alpha=alpha, zorder=1)
476
477
478
479
480
481
482
            self.ax.add_collection(col)
        except ValueError:
            print("Warning: Could not plot ROI circle due to undefined inverse geodesic!")

        self.mark_roi_center()

    def mark_roi_center(self, **kwargs):
483
484
485
486
487
        """
        Mark the ROI center

        :param kwargs: keywords for matplotlib.pyplot.plot() function
        """
488
489
490
491
492
493
494
        kwargs.setdefault('marker', '+')
        kwargs.setdefault('markersize', 20)
        kwargs.setdefault('color', 'k')
        kwargs.setdefault('lw', 2)
        x, y = self.m(self.lon_0, self.lat_0)
        self.m.plot((x), (y), **kwargs)

495
    def plot_grid(self, meridians=None, parallels=None, mer_labels=None, par_labels=None):
496
        """ Plot the longitude and latitude grid in the skypatch """
497
498
499
500
        if meridians is None:
            meridians = np.arange(-180, 181, 60) if abs(self.lat_0) > 60 else np.arange(-180, 181, 20)
        if parallels is None:
            parallels = np.arange(-90, 91, 15) if abs(self.lat_0) > 60 else np.arange(-90, 91, 20)
501

502
503
        self.m.drawmeridians(meridians, labels=[False, False, True, False] if mer_labels is None else mer_labels)
        self.m.drawparallels(parallels, labels=[True, True, False, False] if par_labels is None else par_labels)
504

505
    def plot_thrust(self, n, t, **kwargs):
506
507
508
509
510
511
512
513
        """
        Visualize the thrust observables in the ROI.

        :param n: Thrust axis as given by astrotools.obs.thrust()[1]
        :param t: Thrust values as returned by astrotools.obs.thrust()[0]
        :param kwargs: Keywords passed to matplotlib.pyplot.plot() for axis visualization
        """
        kwargs.setdefault('c', 'red')
514
515
        linestyle_may = kwargs.pop('linestyle', 'solid')
        alpha_may = kwargs.pop('alpha', 0.5)
516
517
518

        lon, lat = coord.vec2ang(n[0])
        # fill thrust array (unit vector phi runs in negative lon direction)
519
520
521
522
        e_phi = coord.sph_unit_vectors(lon, lat)[1]
        sign = np.sign(e_phi[2] - n[1][2])
        phi_major = sign * coord.angle(e_phi, n[1])[0]
        phi_minor = sign * coord.angle(e_phi, n[2])[0]
523
524
525
526
527
528
529
530
        if np.abs(phi_major - phi_minor) < 0.99 * np.pi / 2.:
            phi_minor = 2*np.pi - phi_minor
        t23_ratio = t[1] / t[2]

        # mark the principal axes n3
        u = np.array(np.cos(phi_minor))
        v = -1. * np.array(np.sin(phi_minor))
        urot, vrot, x, y = self.m.rotate_vector(u, v, np.rad2deg(lon), np.rad2deg(lat), returnxy=True)
531
        _phi = np.arctan2(vrot, urot)
532
        s = self.r_roi * (t[1] / 0.15) * self.scale / t23_ratio
533
534
        self.m.plot([x - np.cos(_phi) * s, x + np.cos(_phi)*s],
                    [y - np.sin(_phi) * s, y + np.sin(_phi)*s], linestyle='dashed', alpha=0.5, **kwargs)
535
536
537
538
539

        # mark the principal axes n2
        u = np.array(np.cos(phi_major))
        v = -1. * np.array(np.sin(phi_major))
        urot, vrot, x, y = self.m.rotate_vector(u, v, np.rad2deg(lon), np.rad2deg(lat), returnxy=True)
540
        _phi = np.arctan2(vrot, urot)
541
        s = self.r_roi * (t[1] / 0.15) * self.scale
542
543
        self.m.plot([x - np.cos(_phi) * s, x + np.cos(_phi)*s],
                    [y - np.sin(_phi) * s, y + np.sin(_phi)*s], linestyle=linestyle_may, alpha=alpha_may, **kwargs)
544
545

        # mark the center point
546
        self.m.plot((x), (y), 'o', color=kwargs.pop('c'), markersize=10)
547

548
    def colorbar(self, mappable, cblabel='Energy [eV]', labelsize=12, ticks=None, **kwargs):
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
        """
        Adds a colorbar to a mappable in matplotlib. Replaces matplotlib colorbar() function.
        Use e.g:

        patch = PlotSkyPatch(...)
        mappable = patch.plot_crs(crs)
        patch.colorbar(mappable)

        :param mappable: Mappable in matplotlib.
        :param clabel: Label for the colorbar
        :param ticks: Ticks for the colorbar (either array-like or integer for number of ticks)
        :param kwargs: Keywords passed to matplotlib colorbar() function
        """
        # add a colorbar
        try:
            kwargs.setdefault('location', 'bottom')
            cb = self.m.colorbar(mappable, **kwargs)
            if (ticks is None) or isinstance(ticks, (int, float)):
                vmin, vmax = mappable.get_clim()
                vmin, vmax = smart_round(vmin), smart_round(vmax)
                n_ticks = float(3) if ticks is None else float(ticks)
                step = smart_round((vmax - vmin) / n_ticks, order=1)
                ticks = np.arange(vmin, vmax, step)
                ticks = ticks[(ticks >= vmin) & (ticks <= vmax)]
            cb.set_ticks(ticks)
574
            cb.set_label(cblabel, fontsize=3*labelsize)
575
            t = ['$10^{%.1f}$' % (f) for f in ticks]
576
            cb.ax.set_xticklabels(t, fontsize=int(max(0.8 * 3 * labelsize, 1)))
577
578
579
580
        except KeyError:
            print("Can not plot colorbar on axis.")

    def text(self, x, y, s, **kwargs):
581
        """ Substitudes matplotlib.pyplot.text() function """
582
583
584
585
        kwargs.setdefault('transform', self.ax.transAxes)
        self.ax.text(x, y, s, **kwargs)

    def savefig(self, path, **kwargs):
586
        """ Substitudes matplotlib savefig() function """
587
588
589
        kwargs.setdefault('dpi', 150)
        kwargs.setdefault('bbox_inches', 'tight')
        self.fig.savefig(path, **kwargs)