Reprojecting SunPy maps from one reference frame to another is most commonly done with full-disk images. The built-in implementation of reproject_to() is also capable of reprojecting submaps or cutouts; however, unless the output WCS are carefully adjusted this is just as slow as reprojecting a full-disk image.
To adjust the output WCS to match the size of the reprojected cutout - and therefore to execute reproject without a potentially 4k x 4k array - the CRPIX keywords must be adjusted to account for the new origin in the desired output matrix shape.
This example uses images from November 2021, when Solar Orbiter was close to the Earth and their fields-of-view were quite similar.
import pandas as pd
import sunpy.map
from astropy.coordinates import SkyCoord
from astropy import units as u
from stix_utils import load_SOLO_SPICE
from visible_from_earth import get_observer
from flare_physics_utils import cartesian_diff
peak_utc_corrected=pd.to_datetime('2021-11-02 07:45:11')
load_SOLO_SPICE(peak_utc_corrected)
sobs,swcs=get_observer(peak_utc_corrected,obs='SOLO',sc=Trueout_shape=(1,1))
mcutout=sunpy.map.Map('aia.lev1_euv_12s.2021-11-02T074511Z.171.image._prepped.fits')
mcutout.peek()
To get the Solar Orbiter WCS, it is necessary to first load the proper SPICE kernels. Note: setting out_shape=(1,1) will help later on when determining the extent of the reprojected array shape.
One can then reproject the AIA map to the Solar Orbiter viewpoint using the built-in _reproject_to():
%%timeit
mr=mcutout.reproject_to(swcs)
35.1 s ± 611 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
However, this is very slow! This is because the default size of the Solar Orbiter WCS array matches the full-disk AIA size of (4096,4096). A lot of time is wasted reprojecting NaNs.
To adjust the Solar Orbiter WCS, new CRPIX WCS keywords must be calculated. For this, two quantities are needed: A reference coordinate for the Solar Orbiter coordinate frame, and the new size of the reprojected data array. A pixel scale factor is optional but not yet incorporated.
A reference coordinate is returned by get_observer() if the argument sc=True, otherwise it can be constructed in the usual way of making a SkyCoord at (0”,0”) in the desired reference frame.
Next, find the extent of the new data array after reprojection (without doing the reprojection of course)!
edges_pix = np.concatenate(sunpy.map.map_edges(mcutout))
edges_coord = mcutout.pixel_to_world(edges_pix[:, 0], edges_pix[:, 1])
new_edges_coord = edges_coord.transform_to(sobs)
new_edges_xpix, new_edges_ypix = swcs.world_to_pixel(new_edges_coord)
left, right = np.min(new_edges_xpix), np.max(new_edges_xpix)
bottom, top = np.min(new_edges_ypix), np.max(new_edges_ypix)
If the rotation of either coordinate is not possible, NaN will be returned as either Tx or Ty. In that case, the extent of the new array can be found by transforming all the coordinates (it doesn’t take as long as it sounds) and finding the on-disk boundaries of those coordinates. Note: It may be possible to extend the boundaries to off-disk with no negative consequences to reproject but this hasn’t been tested yet.
def rotated_bounds_on_disk(mcutout,frame):
'''use this to determine extent of rotated shape in case bottom left and top right contain NaN'''
cc=all_coordinates_from_map(mcutout).transform_to(frame)
on_disk=sunpy.map.coordinate_is_on_solar_disk(cc)
on_disk_coordinates=cc[on_disk]
tx = on_disk_coordinates.Tx.value
ty = on_disk_coordinates.Ty.value
return SkyCoord([np.nanmin(tx), np.nanmax(tx)] * u.arcsec,
[np.nanmin(ty), np.nanmax(ty)] * u.arcsec,
frame=smap.coordinate_frame)
if np.isnan([blr.Tx.value,blr.Ty.value]).any() or np.isnan([trr.Tx.value,trr.Ty.value]).any():
blr,trr=rotated_bounds_on_disk(mcutout,sobs.frame)
width, height=cartesian_diff(blr,trr,index=1) #arcsec
out_shape=(int(width.value),int(height.value))
Make a new FITS header/WCS object for input to reproject:
submap_header=sunpy.map.make_fitswcs_header((1,1),sobs)
The output shape (1,1) will be overwritten, as will the translation of the cutout (CRPIX).
submap_header['crpix1'] -= left
submap_header['crpix2'] -= bottom
submap_header['naxis1'] = int(np.ceil(right - left))
submap_header['naxis2'] = int(np.ceil(top - bottom))
rotated_map = mcutout.reproject_to(submap_header)
This results in a speed-up of 12x for this 1000” x 1000” cutout.
%%timeit
rotated_map=smart_reproject(mcutout)
2.74 s ± 116 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
More Examples
Thanks to A. Shih for robust method of determining reprojected array shape.