1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
"""
===========================
Adding an Earth scale image
===========================
This example shows how to plot a map with an image of the Earth added for scale.
"""
import matplotlib.pyplot as plt
from matplotlib.image import AxesImage
from PIL import Image
import astropy.units as u
from astropy.constants import R_earth
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
import sunpy.map
from sunpy.coordinates.utils import solar_angle_equivalency
from sunpy.data import EARTH_IMAGE
from sunpy.data.sample import AIA_193_CUTOUT01_IMAGE
###############################################################################
# We start with a sample AIA image.
cutout_map = sunpy.map.Map(AIA_193_CUTOUT01_IMAGE)
###############################################################################
# We use a (low-resolution) image of Earth that we provide with `sunpy`. You
# can use a different image as desired, and it should have a transparent
# background. We also crop the image tightly so that we can assume that the
# image width/height are equal to the Earth diameter. Finally, we flip the
# image vertically so that it is oriented correctly when plotted with the pixel
# origin in the lower-left corner, which is the convention for maps.
earth = Image.open(EARTH_IMAGE)
earth = earth.crop(earth.getbbox())
earth = earth.transpose(Image.FLIP_TOP_BOTTOM)
##############################################################################
# The first step in plotting the Earth is to convert the Earth's diameter in km
# to arcsec on the plane-of-sky for the time of the observation. We use the
# observer coordinate from the map. For more information about the
# transformation, see :func:`~sunpy.coordinates.utils.solar_angle_equivalency`.
earth_diameter = 2 * R_earth.to(u.arcsec, equivalencies=solar_angle_equivalency(cutout_map.observer_coordinate))
print(f"Earth diameter = {earth_diameter:.2f}")
##############################################################################
# We then create a WCS header so that the Earth will be plotted correctly
# on any underlying map projection. We center the Earth on a specific point in
# the map.
earth_x = 1000 * u.arcsec
earth_y = -200 * u.arcsec
earth_center = SkyCoord(earth_x, earth_y, frame=cutout_map.coordinate_frame)
earth_wcs = sunpy.map.make_fitswcs_header(
(earth.height, earth.width), earth_center,
scale=earth_diameter / ([earth.width, earth.height] * u.pix)
)
##############################################################################
# Now we plot the AIA image and superimpose the Earth for scale. We manually
# create and add the image artist so that the Earth will not be considered
# when autoscaling plot limits. We also add a text label above the Earth.
fig = plt.figure()
ax = fig.add_subplot(projection=cutout_map)
cutout_map.plot(clip_interval=(1, 99.9)*u.percent)
earth_artist = AxesImage(ax, transform=ax.get_transform(WCS(earth_wcs)))
earth_artist.set_data(earth)
ax.add_artist(earth_artist)
ax.text(
earth_x.to_value('deg'), (earth_y + earth_diameter).to_value('deg'),
'Earth to scale', color='white', fontsize=12, horizontalalignment='center',
transform=ax.get_transform('world')
)
plt.show()
|
