Testing Synthetic MOXSI Data#
[158]:
import warnings
from collections import OrderedDict
from scipy.io import readsav
import matplotlib.pyplot as plt
import astropy.wcs
import astropy.units as u
import astropy.constants
from astropy.coordinates import SkyCoord,SpectralCoord
import sunpy
from sunpy.coordinates import get_earth,Helioprojective
from sunpy.image.transform import affine_transform
import numpy as np
import ndcube
import sunpy.map
from astropy.visualization import ImageNormalize,SqrtStretch,LogStretch,quantity_support
from IPython.display import HTML
Define some useful functions.
[27]:
# The following numbers are from Jake and Albert:
CDELT_SPACE = 5.66 * u.arcsec / u.pix
CDELT_WAVE = 55 * u.milliangstrom / u.pix
# The units are from Athiray (and the proposal figure)
BUNIT = 'ph / (pix s)'
[28]:
@u.quantity_input
def construct_rot_matrix(angle:u.deg):
return np.array([[np.cos(angle), -np.sin(angle), np.sin(angle)],
[np.sin(angle), np.cos(angle), -np.cos(angle)],
[0, 0, 1]])
[29]:
def rmatrix_to_pcij(rmatrix):
pcij_keys = {}
for i in range(rmatrix.shape[0]):
for j in range(rmatrix.shape[1]):
pcij_keys[f'PC{i+1}_{j+1}'] = rmatrix[i,j]
return pcij_keys
[30]:
@u.quantity_input
def rotate_image(data, rmatrix, order=4, missing=0.0,):
if order not in range(6):
raise ValueError("Order must be between 0 and 5.")
# Calculate the shape in pixels to contain all of the image data
extent = np.max(np.abs(np.vstack((data.shape @ rmatrix,
data.shape @ rmatrix.T))), axis=0)
# Calculate the needed padding or unpadding
diff = np.asarray(np.ceil((extent - data.shape) / 2), dtype=int).ravel()
# Pad the image array
pad_x = int(np.max((diff[1], 0)))
pad_y = int(np.max((diff[0], 0)))
new_data = np.pad(data,
((pad_y, pad_y), (pad_x, pad_x)),
mode='constant',
constant_values=(missing, missing))
# All of the following pixel calculations use a pixel origin of 0
pixel_array_center = (np.flipud(new_data.shape) - 1) / 2.0
# Apply the rotation to the image data
new_data = affine_transform(new_data.T,
np.asarray(rmatrix),
order=order,
scale=1.0,
image_center=np.flipud(pixel_array_center),
recenter=False,
missing=missing,
use_scipy=True).T
# Unpad the array if necessary
unpad_x = -np.min((diff[1], 0))
if unpad_x > 0:
new_data = new_data[:, unpad_x:-unpad_x]
unpad_y = -np.min((diff[0], 0))
if unpad_y > 0:
new_data = new_data[unpad_y:-unpad_y, :]
return new_data
[31]:
def to_overlappogram(cube, dispersion_angle=0*u.deg, clip=True, order=0):
"""
Flatten intensity cube into an overlappogram such that the latitude and wavelength
directions overlap. The array is clipped such that it starts halfway (in latitude)
through the first wavelength slice and ends halfway through the last wavelength
slice.
TODO: return an NDCube with an appropriate WCS
TODO: generalize such that the dispersion direction can be along any axis or even any
non-orthogonal path through lon,lat space
Parameters
----------
dispersion_angle : `float`, optional
Angle between the dispersion axis and the latitude world coordinate. The dispersion axis
is always aligned with the y-like (row-indexing) pixel axis. A dispersion angle of 0
means that the dispersion is completely latitudinal.
clip : `bool`, optional
If true, ensure that the dispersion direction has the same shape as the wavelength
dimension.
"""
if dispersion_angle % (360*u.deg) == 0:
rot_data = cube.data
else:
rmatrix = construct_rot_matrix(dispersion_angle)
# apply the necessary rotation to every slice in the cube.
# this is an overly complicated way to do the rotation, but
# will suffice for now
layers = []
for d in cube.data:
layers.append(rotate_image(d, rmatrix[:2,:2], order=order))
rot_data = np.array(layers)
shape = rot_data.shape
n_y = int(shape[0] + shape[1])
n_x = int(shape[2])
overlappogram = np.zeros((n_y, n_x))
for i in range(shape[0]):
overlappogram[i:(i+shape[1]), :] += rot_data[i, :, :]
if clip:
# Clip to desired range
# When shape[1] is even, clip_1 == clip_2
# When shape[1] is odd, we arbitrarily shift down by a half index at the short wavelength end
# and up a half at the long wavelength end such that the latitude/wavelength dimension is
# cropped appropriately.
clip_1 = int(np.floor(shape[1]/2))
clip_2 = int(np.ceil(shape[1]/2))
# Add stuff in here to make a WCS???
overlappogram = overlappogram[clip_1:(overlappogram.shape[0]-clip_2),:]
return overlappogram * cube.unit
[32]:
def strided_overlappogram(overlappogram):
"""
Return a "strided" version of the overlappogram array.
For an overlappogram image of shape (N_lam, N_pix), this
function creates an array of dimension (N_lam, N_lam, N_pix)
where each layer is a view of the original array. In other
words, the values at (k,i,j) and (k+1,i,j) point to the
same place in memory.
"""
return np.lib.stride_tricks.as_strided(
overlappogram,
shape=overlappogram.shape[:1]+overlappogram.shape,
strides=(0,)+overlappogram.strides,
writeable=False
)
[33]:
def make_moxsi_ndcube(data, wavelength):
# Here, assume that the data cube has three dimensions and that
# the first dimension corresponds to wavelength, then latitude,
# then longitude.
moxsi_wcs = {
'CRVAL1': 0, # Assume for now that the sun is at the center of the image.
'CRVAL2': 0, # Assume for now that the sun is at the center of the image.
'CRVAL3': wavelength[0],
'CRPIX1': (data.shape[2] + 1) / 2,
'CRPIX2': (data.shape[1] + 1) / 2,
'CRPIX3': 1,
'CDELT1': CDELT_SPACE.to('arcsec / pix').value,
'CDELT2': CDELT_SPACE.to('arcsec / pix').value,
'CDELT3': CDELT_WAVE.to('angstrom / pix').value,
'CUNIT1': 'arcsec',
'CUNIT2': 'arcsec',
'CUNIT3': 'angstrom',
'CTYPE1': 'HPLN-TAN',
'CTYPE2': 'HPLT-TAN',
'CTYPE3': 'WAVE',
}
wcs = astropy.wcs.WCS(moxsi_wcs,)
return ndcube.NDCube(data, wcs=wcs, unit=BUNIT)
Read in some synthetic data
[9]:
#flare_feldman = readsav('../data/forDan_MOXSI_DAATA_X1Class_flare_Feldman.sav')
#flare_scott = readsav('../data/forDan_MOXSI_DAATA_X1Class_flare_scott2015_scaledup.sav')
moxsi_data = readsav('../data/forDan_MOXSI_DATA_09112020_0440_feldman.sav')
#ar_scott = readsav('../data/forDan_MOXSI_DATA_09112020_0440_scott2015_scaledup.sav')
Now make a datacube out of it.
[10]:
moxsi_cube = make_moxsi_ndcube(moxsi_data['moxsi1_img'], moxsi_data['cubixss_wave'])
WARNING: FITSFixedWarning: 'unitfix' made the change 'Changed units:
'angstrom' -> 'Angstrom'. [astropy.wcs.wcs]
[11]:
moxsi_cube[500,...].plot(norm=ImageNormalize(stretch=LogStretch()))
[11]:
<WCSAxesSubplot:>
[12]:
ax = moxsi_cube[:,175,175].plot()
ax.coords[0].set_format_unit(u.angstrom)
ax.coords[0].set_major_formatter('x.x')
Now, the challenge is how to represent this as an overlappogram. We can fairly easily flatten the datacube into something that looks like our overlappogram
[28]:
moxsi_overlappogram_array = to_overlappogram(moxsi_cube,)
[29]:
plt.figure(figsize=(5,15))
plt.imshow(moxsi_overlappogram_array, norm=ImageNormalize(vmax=1e2,stretch=LogStretch()), origin='lower', cmap='hinodexrt')
[29]:
<matplotlib.image.AxesImage at 0x7fc1870e6c10>
We can also plot a “spectra” by summing along the columns.
[30]:
overlap_spectra = moxsi_overlappogram_array.sum(axis=1)
with quantity_support():
plt.plot(moxsi_cube.axis_world_coords(0)[0].to('angstrom'),overlap_spectra)
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
But what about the WCS information? This is just an array with no coordinates. How do we encode this into the WCS?
Below, let’s try a series of pretty bad ideas.
Bad Idea 1: Manually increment the reference pixel#
Our wavelength axis is aligned along our dispersion axis. In this case, the dispersion axis is aligned with latitiude (generally speaking, this will not be true). Thus, when we increment in wavelength, we also want to increment in latitude.
One way to think about this is we want to move the center of the FOV when we slice in wavelength such that the FOV (in lat,lon) is centered on our selected pixel.
In the case we describe above, wavelength only varies with latitude so we can increment the reference pixel corresponding to latitude and leave the longitude pixel untouched.
[31]:
def make_bad_overlappogram_map(cube, k):
wcs = astropy.wcs.WCS({
'CRVAL1': 0,
'CRVAL2': 0,
'CRPIX1': cube.wcs.wcs.crpix[0],
'CRPIX2': k+1,
'CDELT1': cube.wcs.wcs.cdelt[0],
'CDELT2': cube.wcs.wcs.cdelt[1],
'CUNIT1': cube.wcs.wcs.cunit[0].to_string(),
'CUNIT2': cube.wcs.wcs.cunit[1].to_string(),
'CTYPE1': cube.wcs.wcs.ctype[0],
'CTYPE2': cube.wcs.wcs.ctype[1],
})
bad_map = sunpy.map.Map(to_overlappogram(cube), wcs)
bad_map.plot_settings['cmap'] = 'hinodexrt'
bad_map.plot_settings['norm'] = ImageNormalize(stretch=LogStretch())
return bad_map
Let’s choose the first wavelength, \(k=0\)
[32]:
make_bad_overlappogram_map(moxsi_cube,0).peek()
WARNING: SunpyMetadataWarning: Missing metadata for observation time, setting observation time to current time. Set the 'DATE-AVG' FITS keyword to prevent this warning. [sunpy.map.mapbase]
WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
[sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
We can also draw a grid
[33]:
plt.figure(figsize=(3,8))
bm = make_bad_overlappogram_map(moxsi_cube,0)
bm.plot()
_ = bm.draw_grid()
WARNING: SunpyMetadataWarning: Missing metadata for observation time, setting observation time to current time. Set the 'DATE-AVG' FITS keyword to prevent this warning. [sunpy.map.mapbase]
WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
[sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
What if we instead select the wavelength which maximizes our spectra?
[34]:
kmax = np.argmax(overlap_spectra)
plt.figure(figsize=(3,8))
bm = make_bad_overlappogram_map(moxsi_cube,kmax)
bm.plot()
_ = bm.draw_grid()
WARNING: SunpyMetadataWarning: Missing metadata for observation time, setting observation time to current time. Set the 'DATE-AVG' FITS keyword to prevent this warning. [sunpy.map.mapbase]
WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
[sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
Note that we could also crop this image to get only that image that is centered on the wavelength of interest
[35]:
bm_cropped = bm.submap(SkyCoord(-1000,-1000,unit='arcsec',frame=bm.coordinate_frame),
width=2000*u.arcsec,height=2000*u.arcsec)
[36]:
plt.figure(figsize=(8,8))
bm_cropped.plot()
_ = bm_cropped.draw_grid()
WARNING: SunpyMetadataWarning: Missing metadata for observation time, setting observation time to current time. Set the 'DATE-AVG' FITS keyword to prevent this warning. [sunpy.map.mapbase]
WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
[sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
This kind of gives us what we want, but has many problems:
does not work in any general way (what if the dispersion axis is aligned differently?)
separate WCS for each wavelength
No real ability to “slice” in wavelength
Each WCS is separate
Bad Idea 2: Combine these with a NDCubeSequence#
[39]:
def make_slice_wcs(data,wave,k):
return astropy.wcs.WCS({
'CRVAL1': 0,
'CRVAL2': 0,
'CRVAL3': wave[k].to('angstrom').value,
'CRPIX1': (data.shape[2] + 1) / 2,
'CRPIX2': k+1,
'CRPIX3': 1,
'CDELT1': CDELT_SPACE.to('arcsec / pix').value,
'CDELT2': CDELT_SPACE.to('arcsec / pix').value,
'CDELT3': CDELT_WAVE.to('angstrom / pix').value,
'CUNIT1': 'arcsec',
'CUNIT2': 'arcsec',
'CUNIT3': 'angstrom',
'CTYPE1': 'HPLN-TAN',
'CTYPE2': 'HPLT-TAN',
'CTYPE3': 'WAVE',
})
def make_cube_seq(cube):
# Flatten to overlappogram
overlap = to_overlappogram(cube)
# Make strided 3D array
wave = cube.axis_world_coords(0)[0].to('angstrom')
strided_overlap = strided_overlappogram(overlap)
# Construct a "cube" for each wavelength
cubes = []
for i in range(wave.shape[0]):
slice_data = strided_overlap[i:i+1,...]
slice_wcs = make_slice_wcs(slice_data,wave,i)
cubes.append(ndcube.NDCube(slice_data,wcs=slice_wcs,unit=cube.unit))
# Combine them into a sequence where the common axis is wavelength
return ndcube.NDCubeSequence(cubes,common_axis=0)
[40]:
overlap_seq = make_cube_seq(moxsi_cube)
WARNING: FITSFixedWarning: 'unitfix' made the change 'Changed units:
'angstrom' -> 'Angstrom'. [astropy.wcs.wcs]
[41]:
overlap_seq
[41]:
<ndcube.ndcube_sequence.NDCubeSequence object at 0x7fc186d52d00>
NDCubeSequence
--------------
Dimensions: (<Quantity 1073. pix>, <Quantity 1. pix>, <Quantity 1073. pix>, <Quantity 350. pix>)
Physical Types of Axes: [('meta.obs.sequence',), ('em.wl',), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat'), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat')]
Common Cube Axis: 0
[42]:
overlap_seq.cube_like_dimensions
[42]:
$[1073,~1073,~350] \; \mathrm{pix}$
We can now index this object along the first axis such that it returns to us a view of the array, but with the origin shifted according to the selected wavelength.
[43]:
overlap_seq.index_as_cube[kmax,...]
[43]:
<ndcube.ndcube.NDCube object at 0x7fc13d172280>
NDCube
------
Dimensions: [1073. 350.] pix
Physical Types of Axes: [('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat'), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat')]
Unit: ph / (pix s)
Data Type: float64
[44]:
plt.figure(figsize=(3,8))
overlap_seq.index_as_cube[kmax,...].plot(cmap='hinodexrt',norm=ImageNormalize(stretch=LogStretch()))
[44]:
<WCSAxesSubplot:>
[45]:
plt.figure(figsize=(3,8))
overlap_seq.index_as_cube[1073//2,...].plot(cmap='hinodexrt',norm=ImageNormalize(stretch=LogStretch()))
[45]:
<WCSAxesSubplot:>
Note that we haven’t actually duplicated the array \(N_{\lambda}\) times as that would be a massive waste of memory. Instead, we’ve created an array with “0-stride” in the wavelength dimension such that every slice (in wavelength) points to the same place in memory.
There are a few ways to verify this with numpy.
[46]:
p0,_ = overlap_seq.index_as_cube[0].data.__array_interface__['data']
[47]:
p1,_ = overlap_seq.index_as_cube[1].data.__array_interface__['data']
[48]:
p0 == p1
[48]:
True
[49]:
np.shares_memory(overlap_seq.index_as_cube[0].data,
overlap_seq.index_as_cube[1].data)
[49]:
True
[50]:
overlap_seq.index_as_cube[0].data.base is overlap_seq.index_as_cube[1].data.base
[50]:
True
If we can load these into a map sequence, we can also easily animate them (animate NDCubeSequences is a pain).
[51]:
maps = []
wave = moxsi_cube.axis_world_coords(0)[0]
earth = get_earth()
indices = list(range(int(overlap_seq.cube_like_dimensions[0].value)))
for i in indices[::50]:
sl = overlap_seq.index_as_cube[i,...]
m = sunpy.map.Map(sl.data, sl.wcs.low_level_wcs._wcs)
m.meta['WAVELNTH'] = wave[i].to('angstrom').value
m.meta['WAVEUNIT'] = wave.to('angstrom').unit.to_string()
m.meta['DATE-AVG'] = earth.obstime.isot
m.meta['HGLN_OBS'] = earth.lon.to('deg').value
m.meta['HGLT_OBS'] = earth.lat.to('deg').value
m.meta['DSUN_OBS'] = earth.radius.to('m').value
m.meta['RSUN_REF'] = astropy.constants.R_sun.to('m').value
m.plot_settings['cmap'] = 'hinodexrt'
m.plot_settings['norm'] = ImageNormalize(stretch=LogStretch())
maps.append(m)
map_seq = sunpy.map.MapSequence(maps, sortby=None,derotate=False)
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
[52]:
ani = map_seq.plot()
[53]:
HTML(ani.to_jshtml())
[53]:
Slightly Better Idea 3: Encode degeneracy in the WCS#
We can use the PCi_j formalism to couple the wavelength and latitude dimensions. Say we have three pixel dimensions \((p_1,p_2,p_3)\) corresponding to the columns (\(x\)-like), rows (\(y\)-like), and layers (“0-stride” dimension) and three world dimensions \((\theta_x,\theta_y,\lambda)\).
The conversion between the world and pixel coordinates (ignoring for now the complexities of the actual celestial coordinates),
\[\begin{split}\theta_x = \Delta_x(p_1 - r_1),\\
\theta_y = \Delta_y(p_2 - p_3),\\
\lambda = \Delta_\lambda(p_3 - r_3).\\\end{split}\]
Skipping some algebra, this results in the following PCi_j matrix,
\[\begin{split}M =
\begin{bmatrix}
1 & 0 & 0 \\
0 & 1 & -1 \\
0 & 0 & 1
\end{bmatrix}\end{split}\]
This formalism is described in detail in Greisen and Calabretta (2002). Note that because the identity is the default value for PCi_j, we only need specify the off-diagonal elements here, in this case PC2_3=-1.
[55]:
# Flatten to overlappogram
moxsi_overlap = to_overlappogram(moxsi_cube)
# Make strided 3D array
wave = moxsi_cube.axis_world_coords(0)[0].to('angstrom')
moxsi_strided_overlap = strided_overlappogram(moxsi_overlap)
degen_wcs = astropy.wcs.WCS({
'WCSAXES': 3,
'NAXIS1': moxsi_strided_overlap.shape[2],
'NAXIS2': moxsi_strided_overlap.shape[1],
'NAXIS3': moxsi_strided_overlap.shape[0],
'CDELT1': CDELT_SPACE.to('arcsec / pix').value,
'CDELT2': CDELT_SPACE.to('arcsec / pix').value,
'CDELT3': CDELT_WAVE.to('angstrom / pix').value,
'CUNIT1': 'arcsec',
'CUNIT2': 'arcsec',
'CUNIT3': 'Angstrom',
'CTYPE1': 'HPLN-TAN',
'CTYPE2': 'HPLT-TAN',
'CTYPE3': 'WAVE',
'CRPIX1': (moxsi_strided_overlap.shape[2] + 1)/2,
'CRPIX2': (moxsi_strided_overlap.shape[1] + 1)/2,
'CRPIX3': (moxsi_strided_overlap.shape[0] + 1)/2,
'CRVAL1': 0,
'CRVAL2': 0,
'CRVAL3': ((wave[0] + wave[-1])/2).to('angstrom').value,
'PC2_3': -1, # i==j assumed 1, i!=j assumed 0 unless explicitly stated
})
degen_cube = ndcube.NDCube(moxsi_strided_overlap, wcs=degen_wcs)
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
Note that our longitude changes as we move in wavelength.
[56]:
plt.figure(figsize=(3,8))
degen_cube[0,].plot(cmap='hinodexrt',norm=ImageNormalize(stretch=LogStretch()))
[56]:
<WCSAxesSubplot:>
[58]:
plt.figure(figsize=(3,8))
degen_cube[kmax].plot(cmap='hinodexrt',norm=ImageNormalize(stretch=LogStretch()))
[58]:
<WCSAxesSubplot:>
[131]:
plt.figure(figsize=(3,8))
degen_cube[500,].plot(cmap='hinodexrt',norm=ImageNormalize(stretch=LogStretch()))
[131]:
<WCSAxesSubplot:>
Additionally, note that when we plot a slice in wavelength, both the longitude and wavelength axes vary and the alignment between wavelength and latitude varies latitudinally.
[83]:
ax = degen_cube[:,1073//2,175].plot()
ax.coords[2].set_format_unit(u.angstrom)
ax.coords[2].set_major_formatter('x.x')
[57]:
plt.figure(figsize=(3,8))
ani = degen_cube[:50,...].plot(cmap='hinodexrt',norm=ImageNormalize(stretch=LogStretch()))
<Figure size 216x576 with 0 Axes>
[58]:
fani = ani.get_animation()
[60]:
HTML(fani.to_jshtml())
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
[60]:
Possibly Good Idea 4: Generalized PCi_j matrix#
Now the question is, can we come up with a generalized prescription for construction a PCi_j matrix for an arbitrary orientation of the dispersion axis.
[324]:
moxsi_cube.wcs.wcs.cunit
[324]:
['deg', 'deg', 'm']
[328]:
(wave[1]-wave[0]).to('mAngstrom')
[328]:
$55 \; \mathrm{m\mathring{A}}$
[34]:
def construct_generalized_overlappogram(cube, angle):
rmatrix = construct_rot_matrix(angle)
# Flatten to overlappogram
moxsi_overlap = to_overlappogram(cube, dispersion_angle=angle)
# Make strided 3D array
wave = cube.axis_world_coords(0)[0].to('angstrom')
moxsi_strided_overlap = strided_overlappogram(moxsi_overlap)
wcs_keys = {
'WCSAXES': 3,
'NAXIS1': moxsi_strided_overlap.shape[2],
'NAXIS2': moxsi_strided_overlap.shape[1],
'NAXIS3': moxsi_strided_overlap.shape[0],
'CDELT1': CDELT_SPACE.to('arcsec / pix').value,
'CDELT2': CDELT_SPACE.to('arcsec / pix').value,
'CDELT3': CDELT_WAVE.to('angstrom / pix').value,
'CUNIT1': 'arcsec',
'CUNIT2': 'arcsec',
'CUNIT3': 'Angstrom',
'CTYPE1': 'HPLN-TAN',
'CTYPE2': 'HPLT-TAN',
'CTYPE3': 'WAVE',
'CRPIX1': (moxsi_strided_overlap.shape[2] + 1)/2,
'CRPIX2': (moxsi_strided_overlap.shape[1] + 1)/2,
'CRPIX3': (moxsi_strided_overlap.shape[0] + 1)/2,
'CRVAL1': 0,
'CRVAL2': 0,
'CRVAL3': ((wave[0] + wave[-1])/2).to('angstrom').value,
}
pcij_keys = rmatrix_to_pcij(rmatrix)
for k in pcij_keys:
wcs_keys[k] = pcij_keys[k]
wcs = astropy.wcs.WCS(wcs_keys)
overlap_cube = ndcube.NDCube(moxsi_strided_overlap, wcs=wcs)
return overlap_cube
[35]:
moxsi_overlap1 = construct_generalized_overlappogram(moxsi_cube, 0*u.deg)
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
[335]:
plt.figure(figsize=(3,8))
ax = moxsi_overlap1[1073//2].plot(
#plot_axes=('x','y'),
cmap='hinodexrt',
norm=ImageNormalize(stretch=LogStretch())
)
[134]:
ax.coords
[134]:
<CoordinatesMap with 2 world coordinates:
index aliases type unit wrap format_unit visible
----- -------------------------------------------- --------- ---- ----- ----------- -------
0 custom:pos.helioprojective.lon hpln-tan hpln longitude deg 180.0 arcsec yes
1 custom:pos.helioprojective.lat hplt-tan hplt latitude deg None arcsec yes
>
[153]:
moxsi_overlap2 = construct_generalized_overlappogram(moxsi_cube, 90*u.deg)
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
[154]:
moxsi_overlap2.wcs
[154]:
WCS Keywords
Number of WCS axes: 3
CTYPE : 'HPLN-TAN' 'HPLT-TAN' 'WAVE'
CRVAL : 0.0 0.0 3.0480000000000002e-09
CRPIX : 175.5 537.0 537.0
PC1_1 PC1_2 PC1_3 : 6.12323399573676e-17 -1.0 1.0
PC2_1 PC2_2 PC2_3 : 1.0 6.12323399573676e-17 -6.1232339957367e-17
PC3_1 PC3_2 PC3_3 : 0.0 0.0 1.0
CDELT : 0.0015722222222222223 0.0015722222222222223 5.5e-12
NAXIS : 350 1073 1073
[156]:
ax = moxsi_overlap2[10,...].plot(
#plot_axes=('x','y'),
cmap='hinodexrt',
norm=ImageNormalize(stretch=LogStretch())
)
[131]:
ax.coords
[131]:
<CoordinatesMap with 2 world coordinates:
index aliases type unit wrap format_unit visible
----- -------------------------------------------- --------- ---- ----- ----------- -------
0 custom:pos.helioprojective.lon hpln-tan hpln longitude deg 180.0 arcsec yes
1 custom:pos.helioprojective.lat hplt-tan hplt latitude deg None arcsec yes
>
[64]:
moxsi_overlap2[0,1073//2,:].plot()
[64]:
<WCSAxesSubplot:ylabel='Data'>
[47]:
moxsi_overlap3 = construct_generalized_overlappogram(moxsi_cube, 45*u.deg)
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
[61]:
moxsi_overlap3[1073//2].plot(plot_axes=('x','y'),cmap='hinodexrt',norm=ImageNormalize(stretch=LogStretch()))
[61]:
<WCSAxesSubplot:>
[57]:
moxsi_overlap3[0,:,100].plot()
[57]:
<WCSAxesSubplot:ylabel='Data'>
[67]:
moxsi_overlap4 = construct_generalized_overlappogram(moxsi_cube, 5*u.deg)
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
[113]:
plt.figure(figsize=(4,8))
moxsi_overlap4[0].plot(plot_axes=('y','x'),cmap='hinodexrt',norm=ImageNormalize(stretch=LogStretch()))
[113]:
<WCSAxesSubplot:>
[96]:
ax = moxsi_overlap2[:,0,0].plot()
ax.coords[2].set_format_unit(u.angstrom)
ax.coords[2].set_major_formatter('x.x')
[121]:
fig = plt.figure(figsize=(8,8))
foo0 = construct_generalized_overlappogram(moxsi_cube,0*u.deg)
foo = construct_generalized_overlappogram(moxsi_cube, 90*u.deg)
ax = fig.add_subplot(121,projection=foo0[0].wcs)
foo0[0].plot(axes=ax,plot_axes=('y','x'),cmap='hinodexrt',norm=ImageNormalize(stretch=LogStretch()))
ax = fig.add_subplot(122,projection=foo[0].wcs)
foo[0].plot(axes=ax,plot_axes=('y','x'),cmap='hinodexrt',norm=ImageNormalize(stretch=LogStretch()))
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
[121]:
<WCSAxesSubplot:>
[148]:
foo.wcs
[148]:
WCS Keywords
Number of WCS axes: 3
CTYPE : 'HPLN-TAN' 'HPLT-TAN' 'WAVE'
CRVAL : 0.0 0.0 3.0480000000000002e-09
CRPIX : 175.5 537.0 537.0
PC1_1 PC1_2 PC1_3 : 6.12323399573676e-17 -1.0 1.0
PC2_1 PC2_2 PC2_3 : 1.0 6.12323399573676e-17 -6.1232339957367e-17
PC3_1 PC3_2 PC3_3 : 0.0 0.0 1.0
CDELT : 0.0015722222222222223 0.0015722222222222223 5.5e-12
NAXIS : 350 1073 1073
[149]:
foo.wcs.axis_correlation_matrix
[149]:
array([[ True, True, True],
[ True, True, True],
[False, False, True]])
Testing Some Cropping Operations#
[300]:
wave = moxsi_cube.axis_world_coords(0)[0]
iwave_min = 1073//2-5
iwave_max = 1073//2+5
ll = [SkyCoord(-20,-20,unit='arcsec',frame=Helioprojective),SpectralCoord(wave[iwave_min])]
ur = [SkyCoord(20,20,unit='arcsec',frame=Helioprojective),SpectralCoord(wave[iwave_max])]
lr = [SkyCoord(20,-20,unit='arcsec',frame=Helioprojective),SpectralCoord(wave[iwave_max])]
ul = [SkyCoord(-20,20,unit='arcsec',frame=Helioprojective),SpectralCoord(wave[iwave_min])]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
[301]:
moxsi_overlap_cropped = moxsi_overlap1.crop(ll,ur,lr,ul)
WARNING: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change [astropy.wcs.wcsapi.fitswcs]
[302]:
moxsi_overlap_cropped
[302]:
<ndcube.ndcube.NDCube object at 0x7fe7f8319430>
NDCube
------
Dimensions: [11. 19. 8.] pix
Physical Types of Axes: [('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat', 'em.wl'), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat'), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat')]
Unit: None
Data Type: float64
[332]:
moxsi_overlap_cropped[5,...].plot()
[332]:
<WCSAxesSubplot:>
[184]:
moxsi_overlap_cropped[0,:,1].plot()
[184]:
<WCSAxesSubplot:ylabel='Data'>
[289]:
moxsi_overlap1.wcs.world_to_pixel(*ll)
[289]:
(array(170.96643108), array(527.46643107), array(531.))
[298]:
moxsi_overlap1.wcs.world_to_pixel(SkyCoord(0,20,unit='arcsec',frame=Helioprojective),SpectralCoord(wave[iwave_max]))
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change [astropy.wcs.wcsapi.fitswcs]
[298]:
(array(174.5), array(544.53356892), array(541.))
[299]:
moxsi_overlap1[541,537:545,170:179].plot()
[299]:
<WCSAxesSubplot:>
[295]:
#moxsi_overlap1.wcs.world_to_pixel(*ur)
moxsi_overlap1.wcs.world_to_pixel(SkyCoord(20,20,unit='arcsec',frame=Helioprojective),SpectralCoord(wave[iwave_min]))
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change [astropy.wcs.wcsapi.fitswcs]
[295]:
(array(178.03356892), array(534.53356893), array(531.))
[297]:
moxsi_overlap1[531,527:536,170:179].plot()
[297]:
<WCSAxesSubplot:>
[286]:
moxsi_overlap1.wcs.pixel_to_world(*moxsi_overlap1.wcs.world_to_pixel(*ur))
[286]:
[<SkyCoord (Helioprojective: obstime=None, rsun=695700.0 km, observer=None): (Tx, Ty) in arcsec
(20., 20.)>,
<SpectralCoord 3.0755e-09 m>]
[281]:
moxsi_overlap1.wcs.pixel_to_world(178,545,541)
[281]:
[<SkyCoord (Helioprojective: obstime=None, rsun=695700.0 km, observer=None): (Tx, Ty) in arcsec
(19.80999994, 22.6399998)>,
<SpectralCoord 3.0755e-09 m>]
[266]:
for Tx in [0,200,500,800]:
p1,p2,p3 = moxsi_overlap1.wcs.world_to_pixel(
SkyCoord(Tx,0,unit='arcsec',frame=Helioprojective),
SpectralCoord(wave)
)
p1 = p1.astype(int)
p2 = p2.astype(int)
p3 = p3.astype(int)
spec00 = moxsi_overlap1.data[p3,p2,p1]
plt.plot(wave,spec00,label=Tx)
plt.legend()
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change [astropy.wcs.wcsapi.fitswcs]
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change [astropy.wcs.wcsapi.fitswcs]
[266]:
<matplotlib.legend.Legend at 0x7fe7fd8c1d90>
[267]:
p1,p2,p3 = moxsi_overlap1.wcs.world_to_pixel(
SkyCoord(0,0,unit='arcsec',frame=Helioprojective),
SpectralCoord(wave)
)
WARNING: target cannot be converted to ICRS, so will not be set on SpectralCoord [astropy.wcs.wcsapi.fitswcs]
WARNING: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change [astropy.wcs.wcsapi.fitswcs]
[269]:
plt.figure(figsize=(3,8))
ax = moxsi_overlap1[0].plot(
#plot_axes=('x','y'),
cmap='hinodexrt',
norm=ImageNormalize(stretch=LogStretch())
)
[277]:
moxsi_overlap1[0,:,80].plot()
[277]:
<WCSAxesSubplot:ylabel='Data'>
[336]:
m = sunpy.map.Map('/Users/willbarnes/sunpy/data/aia_lev1_131a_2011_02_12t15_00_09_62z_image_lev1.fits')
[337]:
m.wcs.has_celestial
[337]:
True
[338]:
from astropy.wcs.utils import wcs_to_celestial_frame
[339]:
wcs_to_celestial_frame(m.wcs)
[339]:
<Helioprojective Frame (obstime=2011-02-12T15:00:09.620, rsun=696000.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2011-02-12T15:00:09.620, rsun=696000.0 km): (lon, lat, radius) in (deg, deg, m)
(0.00971981, -6.69174749, 1.47661935e+11)>)>
[343]:
wcs_to_celestial_frame(moxsi_overlap1.wcs)
[343]:
<Helioprojective Frame (obstime=None, rsun=695700.0 km, observer=None)>
[344]:
m.wcs.to_header()
[344]:
WCSAXES = 2 / Number of coordinate axes
CRPIX1 = 2047.380005 / Pixel coordinate of reference point
CRPIX2 = 2037.920044 / Pixel coordinate of reference point
PC1_1 = 0.99999706228137 / Coordinate transformation matrix element
PC1_2 = 0.002423928344853 / Coordinate transformation matrix element
PC2_1 = -0.002423928344853 / Coordinate transformation matrix element
PC2_2 = 0.99999706228137 / Coordinate transformation matrix element
CDELT1 = 0.00016686055555556 / [deg] Coordinate increment at reference point
CDELT2 = 0.00016686055555556 / [deg] Coordinate increment at reference point
CUNIT1 = 'deg' / Units of coordinate increment and value
CUNIT2 = 'deg' / Units of coordinate increment and value
CTYPE1 = 'HPLN-TAN' / Coordinate type codegnomonic projection
CTYPE2 = 'HPLT-TAN' / Coordinate type codegnomonic projection
CRVAL1 = 0.0 / [deg] Coordinate value at reference point
CRVAL2 = 0.0 / [deg] Coordinate value at reference point
LONPOLE = 180.0 / [deg] Native longitude of celestial pole
LATPOLE = 0.0 / [deg] Native latitude of celestial pole
MJDREF = 0.0 / [d] MJD of fiducial time
DATE-OBS= '2011-02-12T15:00:09.620' / ISO-8601 time of observation
MJD-OBS = 55604.625111343 / [d] MJD of observation
RSUN_REF= 696000000.0 / [m] Solar radius
DSUN_OBS= 147661935403.95 / [m] Distance from centre of Sun to observer
HGLN_OBS= 0.0097198071673758 / [deg] Stonyhurst heliographic lng of observer
HGLT_OBS= -6.6917474868137 / [deg] Heliographic latitude of observer
[345]:
m.wcs.to_header()
[345]:
WCSAXES = 2 / Number of coordinate axes
CRPIX1 = 2047.380005 / Pixel coordinate of reference point
CRPIX2 = 2037.920044 / Pixel coordinate of reference point
PC1_1 = 0.99999706228137 / Coordinate transformation matrix element
PC1_2 = 0.002423928344853 / Coordinate transformation matrix element
PC2_1 = -0.002423928344853 / Coordinate transformation matrix element
PC2_2 = 0.99999706228137 / Coordinate transformation matrix element
CDELT1 = 0.00016686055555556 / [deg] Coordinate increment at reference point
CDELT2 = 0.00016686055555556 / [deg] Coordinate increment at reference point
CUNIT1 = 'deg' / Units of coordinate increment and value
CUNIT2 = 'deg' / Units of coordinate increment and value
CTYPE1 = 'HPLN-TAN' / Coordinate type codegnomonic projection
CTYPE2 = 'HPLT-TAN' / Coordinate type codegnomonic projection
CRVAL1 = 0.0 / [deg] Coordinate value at reference point
CRVAL2 = 0.0 / [deg] Coordinate value at reference point
LONPOLE = 180.0 / [deg] Native longitude of celestial pole
LATPOLE = 0.0 / [deg] Native latitude of celestial pole
MJDREF = 0.0 / [d] MJD of fiducial time
DATE-OBS= '2011-02-12T15:00:09.620' / ISO-8601 time of observation
MJD-OBS = 55604.625111343 / [d] MJD of observation
RSUN_REF= 696000000.0 / [m] Solar radius
DSUN_OBS= 147661935403.95 / [m] Distance from centre of Sun to observer
HGLN_OBS= 0.0097198071673758 / [deg] Stonyhurst heliographic lng of observer
HGLT_OBS= -6.6917474868137 / [deg] Heliographic latitude of observer
[347]:
foo = sunpy.map.Map(moxsi_cube[1073//2].data,
moxsi_cube[1073//2].wcs.low_level_wcs._wcs)
[351]:
foo.rotate(angle=90*u.deg).peek()
WARNING: SunpyMetadataWarning: Missing metadata for observation time, setting observation time to current time. Set the 'DATE-AVG' FITS keyword to prevent this warning. [sunpy.map.mapbase]
WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
[sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
[353]:
m.rotate(angle=90*u.deg).peek()
WARNING: SunpyUserWarning: Integer input data has been cast to float64. [sunpy.image.transform]
[425]:
import sunpy.data.sample
[426]:
m = sunpy.map.Map(sunpy.data.sample.AIA_094_IMAGE)
[358]:
m.peek(vmin=0,vmax=100)
[380]:
m_rot = m.rotate().rotate(angle=90*u.deg)
m_rot.peek(vmin=0,vmax=1e2)
[417]:
import astropy.units as u
import sunpy.data.sample
import sunpy.map
m = sunpy.map.Map(sunpy.data.sample.AIA_094_IMAGE).rotate()
m_rot = m.rotate(angle=90*u.deg)
fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(121,projection=m)
m.plot(axes=ax,vmin=0,vmax=1e2,annotate=False)
ax = fig.add_subplot(122,projection=m_rot)
m_rot.plot(axes=ax,vmin=0,vmax=1e2,annotate=False)
lon,lat = ax.coords
[418]:
ndcube.NDCube(m.data,m.wcs).plot(norm=ImageNormalize(vmin=0,vmax=1e2,stretch=m.plot_settings['norm'].stretch),
cmap=m.plot_settings['cmap'])
[418]:
<WCSAxesSubplot:>
[419]:
ax = ndcube.NDCube(m_rot.data,m_rot.wcs).plot(
norm=ImageNormalize(vmin=0,vmax=1e2,stretch=m.plot_settings['norm'].stretch),
cmap=m.plot_settings['cmap']
)
#ax.coords[0].set_ticklabel_position('all')
#ax.coords[1].set_ticklabel_position('all')
[420]:
fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(121,projection=m.wcs)
ax.imshow(m.data,norm=ImageNormalize(vmin=0,vmax=1e2,stretch=m.plot_settings['norm'].stretch))
ax = fig.add_subplot(122,projection=m_rot.wcs)
ax.imshow(m_rot.data,norm=ImageNormalize(vmin=0,vmax=1e2,stretch=m.plot_settings['norm'].stretch))
[420]:
<matplotlib.image.AxesImage at 0x7fe8019001c0>
[421]:
fig = plt.figure()
ax = fig.add_subplot(122,projection=m_rot.wcs)
ax.imshow(m_rot.data,norm=ImageNormalize(vmin=0,vmax=1e2,stretch=m.plot_settings['norm'].stretch))
lon,lat = ax.coords
lon.set_ticklabel_position('all')
lon.set_axislabel(ax.get_xlabel(), color='tab:blue')
lon.set_ticklabel('tab:blue')
lat.set_ticklabel_position('all')
lat.set_axislabel(ax.get_ylabel(), color='tab:green')
lat.set_ticklabel('tab:green')
[422]:
m_rot.wcs
[422]:
WCS Keywords
Number of WCS axes: 2
CTYPE : 'HPLN-TAN' 'HPLT-TAN'
CRVAL : 0.00089530541880571 0.00038493926472939
CRPIX : 514.5 514.5
PC1_1 PC1_2 : 6.1444174869866e-17 1.0
PC2_1 PC2_2 : -1.0 6.1185176568363e-17
CDELT : 0.00066744222222222 0.00066744222222222
NAXIS : 1028 1028
[423]:
np.linalg.inv(m_rot.rotation_matrix)
[423]:
array([[ 6.11851766e-17, -1.00000000e+00],
[ 1.00000000e+00, 6.14441749e-17]])
[424]:
m_rot.rotation_matrix
[424]:
array([[ 6.14441749e-17, 1.00000000e+00],
[-1.00000000e+00, 6.11851766e-17]])
[427]:
m_rot = m.rotate(angle=30*u.deg)
[441]:
import sunpy.data.sample
import sunpy.map
m = m = sunpy.map.Map(sunpy.data.sample.AIA_094_IMAGE)
m_rot = m.rotate(angle=30*u.deg)
fig = plt.figure()
ax = fig.add_subplot(111,projection=m)
m.plot(axes=ax,vmin=0,vmax=1e2)
lon,lat = ax.coords
lon.grid(lw=1,color='w',ls='-',alpha=1)
lat.grid(lw=1,color='w',ls='-',alpha=1)
overlay = ax.get_coords_overlay(m_rot.wcs)
overlay[0].grid(lw=1,color='C3',ls='-',alpha=1)
overlay[1].grid(lw=1,color='C3',ls='-',alpha=1)
[ ]: