Commit 0ce1c165 authored by Carl Schaffer's avatar Carl Schaffer
Browse files

Merge branch 'locplot_refactor' into 'master'

Locplot refactor

See merge request !177
parents a74f9526 c14d8fdb
......@@ -37,7 +37,7 @@ def loc_preview(fits_file, fn_out):
time = FitsFile(fits_file).obs_time
with catch_warnings():
warnings.simplefilter("ignore")
fig, ax = make_loc_plot(x.value, y.value, time)
fig, ax = make_loc_plot((x.value, y.value), time)
# plt.tight_layout()
fig.savefig(fn_out, bbox_inches="tight", pad_inches=-0.01, dpi=200)
return fn_out
......
......@@ -20,10 +20,10 @@ import matplotlib.pyplot as plt
import numpy as np
from tqdm import tqdm
from kis_tools.gris.GrisFitsFile import GrisFitsFile
from kis_tools.gris.util import get_observers
from kis_tools.util.locplot import make_loc_plot
from . import GrisFitsFile
from .util import get_observers
from ..util.util import reglob, gris_run_number
from kis_tools.util.util import reglob, gris_run_number
class GrisRun(object):
......@@ -97,7 +97,7 @@ class GrisRun(object):
# not, generate instances
elif isinstance(gris_fits_files[0], str):
gris_fits_files = [
GrisFitsFile.GrisFitsFile(gff) for gff in gris_fits_files
GrisFitsFile(gff) for gff in gris_fits_files
]
# store files sorted by filename
......@@ -468,14 +468,14 @@ class GrisRun(object):
plt.show()
def plot_location(self):
coords = self.calculate_bounding_box()
coords = self.get_fov_corners()
# get median date to better represent the observation
dates = sorted(self.query("obs_time"))
date = dates[len(dates) // 2]
fig, ax = make_loc_plot(
coords[:2], coords[2:], date, uncertainties=self.spatial_uncertainties
coords, date, uncertainties=self.spatial_uncertainties
)
return fig, ax
......@@ -512,6 +512,26 @@ class GrisRun(object):
ymin, ymax = min(y_vals), max(y_vals)
return xmin, xmax, ymin, ymax
def get_fov_corners(self):
"""
Determine the corners of the run's fov. Result is calculated from the first first map in the
observation. If there are subsequent maps, they are ignored.
Returns:
x1,x2,x3,x4 : Coordinate tuples in arcseconds in Helioprojective coordnates
"""
first_map_id = self.maps[0]
map_files = sorted(
list(filter(lambda x: x.mapnumber == first_map_id, self.files)),
key=lambda x: x.filename,
)
first_coords = map_files[0]._coords_from_wcs
last_coords = map_files[-1]._coords_from_wcs
x1 = first_coords[0, :].value
x2 = first_coords[-1, :].value
x3 = last_coords[-1, :].value
x4 = last_coords[0, :].value
return x1, x2, x3, x4
# noinspection PyMethodOverriding
class EmptyGrisRun(GrisRun):
......@@ -536,6 +556,7 @@ class EmptyGrisRun(GrisRun):
if __name__ == "__main__":
# files = glob("/dat/sdc/gris/20150426/level1_split/*_001_???_*.fits")
files = glob("Y:\dat\sdc\gris\20191111\level1_split\*_003_???_*.fits")
files = glob(r"Y:\dat\sdc\gris\20140426\level1_split\*_002_???_*.fits")
gr = GrisRun(files)
fig = gr.plot_location()
props = gr.parse()
......@@ -15,6 +15,7 @@ import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
from matplotlib.ticker import FormatStrFormatter
from ..util.locplot import make_loc_plot
matplotlib.use("Agg")
......@@ -51,7 +52,7 @@ def buildLocPic(filePath):
# check if file has already been created, if so use existing version
if not exists(path_out):
fig, ax = make_loc_plot(x_pos, y_pos, date)
fig, ax = make_loc_plot((x_pos, y_pos), date)
# Save picture
fig.savefig(path_out, bbox_inches="tight", pad_inches=-0.01, dpi=200)
......
......@@ -60,5 +60,8 @@ class TestCoords(unittest.TestCase):
print(msg)
self.assertTrue(coords_match, msg=msg)
# Test that FOV corners are calculated without errors
_ = gr.get_fov_corners()
def test_gris_wcs(self):
wcs = GrisWCSGenerator(gris_structure.raw_files[0])
......@@ -12,23 +12,31 @@ from kis_tools.util.locplot import make_loc_plot
class TestPlotting(unittest.TestCase):
@patch("kis_tools.util.locplot.get_hmi_continuum")
def test_locplot(self, patched_get_hmi):
# Set a dummy image for hmi context
patched_get_hmi.return_value = hmi_context[0]
# Set non-interactive backend, if you want to inspect the plots, turn this off!
matplotlib.use("Agg")
plt.ioff()
date = datetime(2014, 6, 28)
fig, ax = make_loc_plot(0, 0, date)
self.assertEqual(len(ax.patches), 0)
plt.close(fig)
fig, ax = make_loc_plot((-200, 200), (-200, 200), date)
self.assertEqual(len(ax.patches), 1)
plt.close(fig)
fig, ax = make_loc_plot(0, 0, date, uncertainties=(50, 50))
self.assertEqual(len(ax.patches), 3)
plt.close(fig)
fig, ax = make_loc_plot((-200, 200), (-200, 200), date, uncertainties=(50, 50))
self.assertEqual(len(ax.patches), 4)
plt.close(fig)
test_coords = [
(0, 0), # X,Y
((-200, 200), (200, -200)), # Bounding Box
((-200, 0), (0, 200), (200, 0), (0, -200)) # FOV corners
]
# Length of patches in axes is used to check for correct plotting, these are the expected
# lengths for the different coordinates
target_lengths = [0, 1, 1]
# loop over inputs and perform plotting
for coords, target_length in zip(test_coords, target_lengths):
# Test plotting
fig, ax = make_loc_plot(coords, date)
self.assertEqual(len(ax.patches), target_length)
plt.close(fig)
# Test plotting with uncertainties
fig, ax = make_loc_plot(coords, date, uncertainties=(50, 50))
self.assertEqual(len(ax.patches), target_length + 3)
plt.close(fig)
......@@ -8,7 +8,6 @@ Created by schaffer at 4/17/19
import datetime
import os
from collections import Iterable
from os.path import exists
from os.path import join
from shutil import move
......@@ -17,7 +16,7 @@ import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import sunpy.map
from matplotlib.patches import Ellipse, Rectangle
from matplotlib.patches import Ellipse, Polygon
from sunpy.net.helioviewer import HelioviewerClient
......@@ -45,94 +44,100 @@ def get_hmi_continuum(date):
return path
def make_loc_plot(X, Y, date, uncertainties=None):
def make_loc_plot(coords, date, uncertainties=None):
"""
Overplot coordionates onto an HMI image.
Args:
X (float or tuple (X_min, X_max)): helioprojective x-coordinate in arcsec
Y (float or tuple (Y_min, Y_max)): helioprojective y-coordinate in arcsec
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)
date (datetime.datetime): date of observation
uncertainties: (sX, sY): coordinate certainty (standard deviation) in arcseconds, defaults to None
Returns:
"""
add_extent = False
if isinstance(X, Iterable) and isinstance(X, Iterable):
add_extent = True
width = max(X) - min(X)
height = max(Y) - min(Y)
X = sum(X) / len(X)
Y = sum(Y) / len(Y)
r = Rectangle(
(X - width / 2, Y - height / 2),
width,
height,
edgecolor="red",
facecolor="none",
)
if str(X) == "-0.0":
X = X * -1
if str(Y) == "-0.0":
Y = Y * -1
# Get SunPicture from Helioviewer.org
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
continuum_file = get_hmi_continuum(date)
# Set position of the labels depending on the position
if X > 0:
xlabelPos = X - 275
ylabelPos = Y - 19
if x_center > 0:
xlabel_pos = x_center - 275
ylabel_pos = y_center - 19
if X < 0:
xlabelPos = X + 275
ylabelPos = Y - 19
if x_center < 0:
xlabel_pos = x_center + 275
ylabel_pos = y_center - 19
if Y > 0:
ylabelPos = Y - 110
xlabelPos = X
if y_center > 0:
ylabel_pos = y_center - 110
xlabel_pos = x_center
if Y < 0:
ylabelPos = Y + 80
xlabelPos = X
if y_center < 0:
ylabel_pos = y_center + 80
xlabel_pos = x_center
if Y == 0 and X == 0:
ylabelPos = Y - 110
xlabelPos = X
if y_center == 0 and x_center == 0:
ylabel_pos = y_center - 110
xlabel_pos = x_center
else:
pass
# Convert position from m in arcsec
xc = [X] * u.arcsec
yc = [Y] * u.arcsec
xc = [x_center] * u.arcsec
yc = [y_center] * u.arcsec
# Create the picture
smap = sunpy.map.Map(continuum_file)
fig, ax = plt.subplots()
smap._default_heliographic_longitude = 0
smap.plot(ax)
hmi_map = sunpy.map.Map(continuum_file)
figure, axes = plt.subplots()
hmi_map._default_heliographic_longitude = 0
hmi_map.plot(axes)
# Add labels
title = ax.title.get_text()
ax.annotate(title, xy=(0, 1010), ha="center", color="white", size=7)
plt.plot(xc, yc, color="red", marker="x", markersize=10, mew=2)
title = axes.title.get_text()
axes.annotate(title, xy=(0, 1010), ha="center", color="white", size=7)
# Plot marker if no corners are given
if corners is None:
plt.plot(xc, yc, color="red", marker="x", markersize=10, mew=2)
plt.tight_layout(pad=0)
marker_label = f"{X:-7.1f}|{Y:-7.1f}"
ax.annotate(
marker_label = f"{x_center:-7.1f}|{y_center:-7.1f}"
axes.annotate(
marker_label,
xy=(xlabelPos, ylabelPos),
xy=(xlabel_pos, ylabel_pos),
ha="center",
color="red",
size=10,
)
# Disable automatic labels
ax.set_title("")
axes.set_title("")
plt.xlabel("")
plt.ylabel("")
# Define x & y axis
ax.set_xlim([-1050, 1050])
ax.set_ylim([-1050, 1050])
ax.tick_params(
axes.set_xlim([-1050, 1050])
axes.set_ylim([-1050, 1050])
axes.tick_params(
axis="y",
direction="in",
pad=-1,
......@@ -143,7 +148,7 @@ def make_loc_plot(X, Y, date, uncertainties=None):
grid_linewidth=0.5,
grid_color="red",
)
ax.tick_params(
axes.tick_params(
axis="x",
direction="in",
pad=-4.5,
......@@ -156,24 +161,24 @@ def make_loc_plot(X, Y, date, uncertainties=None):
)
plt.xticks(np.arange(-1000, 1001, 200))
plt.yticks(np.arange(-1000, 1001, 200))
for tick in ax.yaxis.get_major_ticks():
for tick in axes.yaxis.get_major_ticks():
tick.label1.set_horizontalalignment("left")
tick.label1.set_verticalalignment("bottom")
for tick in ax.xaxis.get_major_ticks():
for tick in axes.xaxis.get_major_ticks():
tick.label1.set_horizontalalignment("left")
tick.label1.set_verticalalignment("bottom")
plt.grid(True)
edge_length = 8
fig.set_size_inches(edge_length, edge_length)
ax.set_facecolor("black")
fig.patch.set_facecolor("black")
figure.set_size_inches(edge_length, edge_length)
axes.set_facecolor("black")
figure.patch.set_facecolor("black")
# Add uncertainty marker if uncertainties are given.
if uncertainties:
one_sigma = Ellipse(
(X, Y),
(x_center, y_center),
2 * uncertainties[0],
2 * uncertainties[1],
edgecolor="red",
......@@ -182,7 +187,7 @@ def make_loc_plot(X, Y, date, uncertainties=None):
label="1 Std. Deviation",
)
two_sigma = Ellipse(
(X, Y),
(x_center, y_center),
2 * 2 * uncertainties[0],
2 * 2 * uncertainties[1],
edgecolor="red",
......@@ -191,7 +196,7 @@ def make_loc_plot(X, Y, date, uncertainties=None):
label="2 Std. Deviations",
)
three_sigma = Ellipse(
(X, Y),
(x_center, y_center),
3 * 2 * uncertainties[0],
3 * 2 * uncertainties[1],
edgecolor="red",
......@@ -199,16 +204,17 @@ def make_loc_plot(X, Y, date, uncertainties=None):
alpha=0.033,
label="3 Std. Deviations",
)
ax.add_patch(one_sigma)
ax.add_patch(two_sigma)
ax.add_patch(three_sigma)
axes.add_patch(one_sigma)
axes.add_patch(two_sigma)
axes.add_patch(three_sigma)
# Add range box if range was specified
if add_extent:
ax.add_patch(r)
return fig, ax
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
if __name__ == "__main__":
fig, ax = make_loc_plot(200, -152, datetime.datetime.now())
fig, ax = make_loc_plot((200, -152), datetime.datetime.now())
plt.show()
......@@ -9,7 +9,7 @@ with open("README.md") as f:
README = f.read()
info = get_git_info(force=True)
version = "3.3.3"
version = "3.4.0"
with open('requirements.txt', 'r') as req_file:
requirements = [line.strip() for line in req_file]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment