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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
|
"""
==============================================================================
Create a Helioprojective Map from observations in the RA-DEC coordinate system
==============================================================================
How to create a `~sunpy.map.Map` in Helioprojective Coordinate Frame from radio observations
in GCRS (RA-DEC).
In this example a LOFAR FITS file (created with LOFAR's `Default Pre-Processing Pipeline (DPPP) and
WSClean Imager <https://support.astron.nl/LOFARImagingCookbook/dppp.html>`__) is read in,
the WCS header information is then used to make a new header with the information in Helioprojective,
and a `~sunpy.map.Map` is made.
The LOFAR example file has a WCS in celestial coordinates i.e. Right Ascension and
Declination (RA-DEC). For this example, we are assuming that the definition of LOFAR's
coordinate system for this observation is exactly the same as Astropy's ~astropy.coordinates.GCRS.
For many solar studies we may want to plot this data in some Sun-centered coordinate frame,
such as `~sunpy.coordinates.frames.Helioprojective`. In this example we read the data and
header information from the LOFAR FITS file and then create a new header with updated WCS
information to create a `~sunpy.map.Map` with a HPC coordinate frame. We will make use of the
`astropy.coordinates` and `sunpy.coordinates` submodules together with
:func:`~sunpy.map.header_helper.make_fitswcs_header` to create a new header and generate a
`~sunpy.map.Map`.
"""
import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
from astropy.coordinates import EarthLocation, SkyCoord
from astropy.io import fits
from astropy.time import Time
import sunpy.data.sample
import sunpy.map
from sunpy.coordinates import frames, sun
##############################################################################
# We will first begin be reading in the header and data from the FITS file.
hdu = fits.open(sunpy.data.sample.LOFAR_IMAGE)
header = hdu[0].header
#####################################################################################
# The data in this file is in a datacube structure to hold difference frequencies and
# polarizations. We are only interested in the image at one frequency (the only data
# in the file) so we index the data array to be the 2D data of interest.
data = hdu[0].data[0, 0, :, :]
################################################################################
# We can inspect the header, for example, we can print the coordinate system and
# type projection, which here is RA-DEC.
print(header['ctype1'], header['ctype2'])
###############################################################################
# Lets pull out the observation time and wavelength from the header, we will use
# these to create our new header.
obstime = Time(header['date-obs'])
frequency = header['crval3']*u.Hz
###############################################################################
# To create a new `~sunpy.map.Map` header we need convert the reference coordinate
# in RA-DEC (that is in the header) to Helioprojective. To do this we will first create
# an `astropy.coordinates.SkyCoord` of the reference coordinate from the header information.
# We will need the location of the observer (i.e. where the observation was taken).
# We first establish the location on Earth from which the observation takes place, in this
# case LOFAR observations are taken from Exloo in the Netherlands, which we define in lat and lon.
# We can convert this to a SkyCoord in GCRSat the observation time.
lofar_loc = EarthLocation(lat=52.905329712*u.deg, lon=6.867996528*u.deg)
lofar_gcrs = SkyCoord(lofar_loc.get_gcrs(obstime))
##########################################################################
# We can then define the reference coordinate in terms of RA-DEC from the header information.
# Here we are using the ``obsgeoloc`` keyword argument to take into account that the observer is not
# at the center of the Earth (i.e. the GCRS origin). The distance here is the Sun-observer distance.
reference_coord = SkyCoord(header['crval1']*u.Unit(header['cunit1']),
header['crval2']*u.Unit(header['cunit2']),
frame='gcrs',
obstime=obstime,
obsgeoloc=lofar_gcrs.cartesian,
obsgeovel=lofar_gcrs.velocity.to_cartesian(),
distance=lofar_gcrs.hcrs.distance)
##########################################################################
# Now we can convert the ``reference_coord`` to the HPC coordinate frame.
reference_coord_arcsec = reference_coord.transform_to(frames.Helioprojective(observer=lofar_gcrs))
##########################################################################
# Now we need to get the other parameters from the header that will be used
# to create the new header - here we can get the cdelt1 and cdelt2 which are
# the spatial scales of the data axes.
cdelt1 = (np.abs(header['cdelt1'])*u.deg).to(u.arcsec)
cdelt2 = (np.abs(header['cdelt2'])*u.deg).to(u.arcsec)
##################################################################################
# Finally, we need to specify the orientation of the HPC coordinate grid because
# GCRS north is not in the same direction as HPC north. For convenience, we use
# :func:`~sunpy.coordinates.sun.P` to calculate this relative rotation angle,
# although due to subtleties in definitions, the returned value is inaccurate by
# 2 arcmin, equivalent to a worst-case shift of 0.6 arcsec for HPC coordinates
# on the disk. The image will need to be rotated by this angle so that solar north
# is pointing up.
P1 = sun.P(obstime)
##########################################################################
# Now we can use this information to create a new header using the helper
# function :func:`~sunpy.map.header_helper.make_fitswcs_header`. This will
# create a MetaDict which we contain all the necessary WCS information to
# create a `~sunpy.map.Map`. We provide a reference coordinate (in HPC),
# the spatial scale of the observation (i.e., ``cdelt1`` and ``cdelt2``),
# and the rotation angle (P1). Note that here, 1 is subtracted from the
# ``crpix1`` and ``crpix2`` values, this is because the ``reference_pixel``
# keyword in :func:~sunpy.map.header_helper.make_fitswcs_header` is zero
# indexed rather than the fits convention of 1 indexed.
new_header = sunpy.map.make_fitswcs_header(data, reference_coord_arcsec,
reference_pixel=u.Quantity([header['crpix1']-1,
header['crpix2']-1]*u.pixel),
scale=u.Quantity([cdelt1, cdelt2]*u.arcsec/u.pix),
rotation_angle=-P1,
wavelength=frequency.to(u.MHz).round(2),
observatory='LOFAR')
##########################################################################
# Let's inspect the new header.
print(new_header)
##########################################################################
# Lets create a `~sunpy.map.Map`.
lofar_map = sunpy.map.Map(data, new_header)
##########################################################################
# We can now plot this map and inspect it.
fig = plt.figure()
ax = fig.add_subplot(projection=lofar_map)
lofar_map.plot(axes=ax, cmap='viridis')
lofar_map.draw_limb(axes=ax)
plt.show()
##########################################################################
# We can now rotate the image so that solar north is pointing up and create
# a submap in the field of view of interest.
lofar_map_rotate = lofar_map.rotate()
bl = SkyCoord(-1500*u.arcsec, -1500*u.arcsec, frame=lofar_map_rotate.coordinate_frame)
tr = SkyCoord(1500*u.arcsec, 1500*u.arcsec, frame=lofar_map_rotate.coordinate_frame)
lofar_submap = lofar_map_rotate.submap(bl, top_right=tr)
##########################################################################
# Now lets plot this map, and overplot some contours.
fig = plt.figure()
ax = fig.add_subplot(projection=lofar_submap)
lofar_submap.plot(axes=ax, cmap='viridis')
lofar_submap.draw_limb(axes=ax)
lofar_submap.draw_grid(axes=ax)
lofar_submap.draw_contours(np.arange(30, 100, 5)*u.percent, axes=ax)
plt.show()
|
