The actual rotation in Map.rotate() is performed by the affine_transform() function in sunpy.image.transform. This function takes a numpy array and rotates it based on a given rotation matrix. It can also scale and shift the array.
Map.rotate() defines the neccessary parameters and updates the appropriate Map header values after the rotation. This does most of the work required for aiaprep(), which simply defines a scale factor and calls Map.rotate().
Each of these functions will be demonstrated below.
affine_transform() rotates the given image using skimage.transform.AffineTransform and skimage.transform.warp(). The rotation matrix corresponding to the desired transformation is given as input to affine_transform()
In [ ]:
# Import statements etc. Nothing unusual here.
import numpy as np
import sunpy
from sunpy.image.transform import affine_transform
from sunpy.instr.aia import aiaprep
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from skimage.data import camera
%matplotlib inline
# Open demo image.
testdata = camera().astype('float')
# Rotate the image 90 degrees anti-clockwise around the centre of the array.
rotdata = affine_transform(testdata, np.array([[0, -1], [1, 0]]))
# Plot the original map against the rotated image.
fig = plt.figure(figsize=(12, 9))
fig.add_subplot(1, 2, 1)
plt.imshow(testdata, cmap=cm.Greys_r)
plt.title('Original image')
fig.add_subplot(1, 2, 2)
plt.imshow(rotdata, cmap=cm.Greys_r)
plt.title('Rotated using affine_transform()')
plt.show()
It can also rotate the image around a user-defined pixel position (which can be moved to the centre of the image or left where it is) and scale the image.
In [ ]:
# Rotate the image by 90 degrees, scale it to twice the original size and move the user-defined
# image centre to the centre of the array.
rotdata = affine_transform(testdata, np.array([[0, -1], [1, 0]]), scale=2.0,
image_center=(400, 256), recenter=True)
# Plot the original map against the rotated image.
fig = plt.figure(figsize=(12, 9))
fig.add_subplot(1, 2, 1)
plt.imshow(testdata, cmap=cm.Greys_r)
plt.title('Original image')
fig.add_subplot(1, 2, 2)
plt.imshow(rotdata, cmap=cm.Greys_r)
plt.title('Rotated using affine_transform()')
plt.show()
Map.rotate() can be given either an angle or a rotation matrix (with the angle and rmatrix keywords, respectively), which will be passed on to affine_transform() (an angle will first be converted into the corresponding rotation matrix). When called without arguments, the map data is rotated using the map's rotation_matrix attribute (which can be calculated using the 'crota2' keyword if the attribute doesn't already exist). This will align the map so that solar North is at the top of the image.
The scale, image_center and recenter keywords can all be used with Map.rotate() to pass those values to affine_transform().
As mentioned above, Map.rotate() can be called with no arguments. In this case, the default behaviour is just to account for the rotation of the instrument.
In [ ]:
# Open image from AIA during an SDO roll maneuver.
aiamap = sunpy.map.Map('/home/stuart/Dropbox/Stu and Lenny/aiacalibim5.fits')
# Rotate your owl.
rotmap = aiamap.rotate()
# Plot original image against the rotated image
fig = plt.figure(figsize=(12, 9))
fig.add_subplot(1, 2, 1)
aiamap.plot()
plt.title('Original image')
fig.add_subplot(1, 2, 2)
rotmap.plot()
plt.title('Image corrected for AIA rotation')
plt.show()
The original map used in the example above has a rotation angle of ~-45 degrees. The corresponding rotation matrix is stored in Map.rotation_matrix.
In [ ]:
# Print rotation info.
print aiamap.rotation_matrix
print aiamap.meta['crota2']
This information is automatically updated by Map.rotate() to account for the new orientation of the image.
In [ ]:
# Print things again.
print rotmap.rotation_matrix
print rotmap.meta['crota2']
Note that the 'crota2' keyword has been removed from the map metadata after rotation.
Although it is being demonstrated here in the context of its uses for aiaprep(), Map.rotate() can be used with any Map with WCS header data. For example, STEREO data:
In [ ]:
# Open image from EUVI during a STEREO roll maneuver.
stereomap = sunpy.map.Map('/home/stuart/Dropbox/Stu and Lenny/stereo_rot_demo.fts')
# Rotate the image by 90 degrees around the centre of the array.
rotmap = stereomap.rotate()
# Plot the original map against the rotated image.
fig = plt.figure(figsize=(12, 9))
fig.add_subplot(1, 2, 1)
stereomap.plot(norm=aiamap.mpl_color_normalizer)
plt.title('Original image')
fig.add_subplot(1, 2, 2)
rotmap.plot(norm=aiamap.mpl_color_normalizer)
plt.title('Image corrected for STEREO rotation')
plt.show()
aiaprep() is the SunPy implementation of IDL's aia_prep() function, which processes an AIA image from level 1 to level 1.5. This processing consists of rotating the image to align solar North with the top of the image, moving the solar centre to the centre of the image, and rescaling the image so that each pixel is exactly 0.6 arcsec square.
This function WILL NOT give the same results as the IDL version, and it isn't intended to - the functions serve the same purpose but use different transformations to process the image.
aiaprep() calls Map.rotate() with the default rotation to align solar North to the top of the image. It also moves solar centre to the centre of the image using the 'recenter' keyword and calculates the factor required to rescale the image to 0.6 arcsec/pixel. The correction for rotation can be seen by using the same AIA image as in the example as above:
In [ ]:
# Process the AIA image
prepmap = aiaprep(aiamap)
# Plot original image against the processed image
fig = plt.figure(figsize=(12, 9))
fig.add_subplot(1, 2, 1)
aiamap.plot()
plt.title('Original image')
fig.add_subplot(1, 2, 2)
prepmap.plot()
plt.title('Processed image')
plt.show()
Similarly, the correction for the instrument offset can be demonstrated by looking at an image taken during a cross maneuver:
In [ ]:
# Open image from AIA during an SDO cross maneuver.
aiamap = sunpy.map.Map('/home/stuart/Dropbox/Stu and Lenny/aiacalibim1.fits')
# Process the AIA image
prepmap = aiaprep(aiamap)
# Plot original image against the processed image
fig = plt.figure(figsize=(12, 9))
fig.add_subplot(1, 2, 1)
aiamap.plot()
aiamap.draw_grid()
plt.title('Original image')
fig.add_subplot(1, 2, 2)
prepmap.plot()
prepmap.draw_grid()
plt.title('Processed image')
plt.show()
Notice that the latitude/longitude grid has been correctly drawn on the processed image, demonstrating that the WCS header information has been updated to account for the shift of the image.