locplot.py 6.07 KB
Newer Older
1
2
3
4
5
6
7
8
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
Created by schaffer at 4/17/19

"""

9
import datetime
10
import os
11
12
13
14
from os.path import exists
from os.path import join
from shutil import move

15
import astropy.units as u
16
import matplotlib.pyplot as plt
17
import numpy as np
18
import sunpy.map
19
from matplotlib.patches import Ellipse, Polygon
20
21
22
23
24
25
26
27
28
29
30
31
32
33
from sunpy.net.helioviewer import HelioviewerClient


def get_hmi_continuum(date):
    """
    Get HMI continuum jp2 image. Check /dat/sdc/hmi_context_cache first, if nothing is found, download from SDO via Helioviewer
    Args:
        date: datetime query

    Returns:
        path: path to the jp2 file
    """
    cache = "/dat/sdc/hmi_context_cache/"
    filename = f"HMI_continuum_{date.strftime('%Y%m%d_%H-%M-%S')}.jp2"
34
35
    if not exists(cache):
        os.makedirs(cache)
36
37
38
39
    path = join(cache, filename)
    if not exists(path):
        hv = HelioviewerClient()
        file = hv.download_jp2(date, source_id=18)
40
41
        if not file:
            raise IOError(f"Could not download file for {date}")
42
        move(file, path)
43

44
    return path
45
46


47
def make_loc_plot(coords, date, uncertainties=None):
48
    """
49
    Overplot coordionates onto an HMI image.
50
51

    Args:
52
53
54
55
        coords: helioprojective coordinates in arcsec, possible shapes are:
                (2) (x_center/y_center)          : Simple point
                (2,2) (point, x_center/y_center) : Bounding Box (upper left and lower right corner)
                (4,2) (point, x_center/y_center) : FOV (Corners of field of view)
56
        date (datetime.datetime):  date of observation
57
        uncertainties: (sX, sY): coordinate certainty (standard deviation) in arcseconds, defaults to None
58
59
60
61
    Returns:


    """
62

63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
    global x_center
    corners = None
    coords = np.array(coords)

    if coords.shape == (2,):
        x_center, y_center = coords
    elif coords.shape == (2, 2):
        x_center, y_center = coords.mean(axis=0)
        xmin, ymin = coords.min(axis=1)
        xmax, ymax = coords.max(axis=1)
        corners = np.array(((xmin, ymin), (xmin, ymax), (xmax, ymax), (xmax, ymin)))
    elif coords.shape == (4, 2):
        x_center, y_center = coords.mean(axis=0)
        corners = coords
    else:
        raise ValueError('Invalid input shape! Use (2,), or (2,2)/(4,2) (Point,X/Y)')

    # Get hmi reference image
81
    continuum_file = get_hmi_continuum(date)
82

83
    # Set position of the labels depending on the position
84
85
86
    if x_center > 0:
        xlabel_pos = x_center - 275
        ylabel_pos = y_center - 19
87

88
89
90
    if x_center < 0:
        xlabel_pos = x_center + 275
        ylabel_pos = y_center - 19
91

92
93
94
    if y_center > 0:
        ylabel_pos = y_center - 110
        xlabel_pos = x_center
95

96
97
98
    if y_center < 0:
        ylabel_pos = y_center + 80
        xlabel_pos = x_center
99

100
101
102
103
104
    if y_center == 0 and x_center == 0:
        ylabel_pos = y_center - 110
        xlabel_pos = x_center
    else:
        pass
105
106

    # Convert position from m in arcsec
107
108
    xc = [x_center] * u.arcsec
    yc = [y_center] * u.arcsec
109
110

    # Create the picture
111
112
113
114
    hmi_map = sunpy.map.Map(continuum_file)
    figure, axes = plt.subplots()
    hmi_map._default_heliographic_longitude = 0
    hmi_map.plot(axes)
115
116

    # Add labels
117
118
    title = axes.title.get_text()
    axes.annotate(title, xy=(0, 1010), ha="center", color="white", size=7)
119
120
121
    # Plot marker if no corners are given
    if corners is None:
        plt.plot(xc, yc, color="red", marker="x", markersize=10, mew=2)
122
    plt.tight_layout(pad=0)
123
124
    marker_label = f"{x_center:-7.1f}|{y_center:-7.1f}"
    axes.annotate(
125
        marker_label,
126
        xy=(xlabel_pos, ylabel_pos),
127
128
129
130
131
132
        ha="center",
        color="red",
        size=10,
    )

    # Disable automatic labels
133
    axes.set_title("")
134
135
136
137
    plt.xlabel("")
    plt.ylabel("")

    # Define x & y axis
Carl Schaffer's avatar
Carl Schaffer committed
138
139
    axes.set_xlim([-1100, 1100])
    axes.set_ylim([-1100, 1100])
140
    axes.tick_params(
141
142
        axis="y",
        direction="in",
Carl Schaffer's avatar
Carl Schaffer committed
143
        pad=-5,
144
        length=0,
Carl Schaffer's avatar
Carl Schaffer committed
145
        labelsize=6,
146
147
148
149
150
        colors="red",
        grid_linestyle="dotted",
        grid_linewidth=0.5,
        grid_color="red",
    )
151
    axes.tick_params(
152
153
        axis="x",
        direction="in",
Carl Schaffer's avatar
Carl Schaffer committed
154
        pad=-5,
155
        length=0,
Carl Schaffer's avatar
Carl Schaffer committed
156
        labelsize=6,
157
158
159
160
161
162
163
        colors="red",
        grid_linestyle="dotted",
        grid_linewidth=0.5,
        grid_color="red",
    )
    plt.xticks(np.arange(-1000, 1001, 200))
    plt.yticks(np.arange(-1000, 1001, 200))
164
    for tick in axes.yaxis.get_major_ticks():
165
166
        tick.label1.set_horizontalalignment("left")
        tick.label1.set_verticalalignment("bottom")
167
    for tick in axes.xaxis.get_major_ticks():
168
169
170
171
        tick.label1.set_horizontalalignment("left")
        tick.label1.set_verticalalignment("bottom")
    plt.grid(True)

172
    edge_length = 8
173

174
175
176
    figure.set_size_inches(edge_length, edge_length)
    axes.set_facecolor("black")
    figure.patch.set_facecolor("black")
177

178
179
    # Add uncertainty marker if uncertainties are given.
    if uncertainties:
Carl Schaffer's avatar
Carl Schaffer committed
180
        one_sigma = Ellipse(
181
            (x_center, y_center),
Carl Schaffer's avatar
Carl Schaffer committed
182
183
184
185
186
187
188
189
            2 * uncertainties[0],
            2 * uncertainties[1],
            edgecolor="red",
            facecolor="red",
            alpha=0.1,
            label="1 Std. Deviation",
        )
        two_sigma = Ellipse(
190
            (x_center, y_center),
Carl Schaffer's avatar
Carl Schaffer committed
191
192
193
194
195
196
197
198
            2 * 2 * uncertainties[0],
            2 * 2 * uncertainties[1],
            edgecolor="red",
            facecolor="red",
            alpha=0.066,
            label="2 Std. Deviations",
        )
        three_sigma = Ellipse(
199
            (x_center, y_center),
Carl Schaffer's avatar
Carl Schaffer committed
200
201
202
203
204
205
206
            3 * 2 * uncertainties[0],
            3 * 2 * uncertainties[1],
            edgecolor="red",
            facecolor="red",
            alpha=0.033,
            label="3 Std. Deviations",
        )
207
208
209
        axes.add_patch(one_sigma)
        axes.add_patch(two_sigma)
        axes.add_patch(three_sigma)
210

211
    # Add range box if range was specified
212
213
214
215
    if corners is not None:
        pol = Polygon(corners, closed=True, label='Field of View', fill=False, color='red')
        axes.add_patch(pol)
    return figure, axes
216
217
218


if __name__ == "__main__":
219
    fig, ax = make_loc_plot((200, -152), datetime.datetime.now())
220
    plt.show()