problem while overplotting limb flares

106 views
Skip to first unread message

Sudha

unread,
Mar 22, 2022, 5:23:57 AM3/22/22
to SunPy
Hi,

I'm using the following versions of sunpy and astropy
print(sunpy.__version__)
3.1.5
print(astropy.__version__)
5.0.2
print(matplotlib.__version__)
3.5.1


I have SDO AIA 193 and Hinode XRT Be Thick fits files for a limb event. My aim is to over-plot the two and I seem to have hit a problem.

The steps I've followed are as follows:

Tested composite map creation and over-plotting using the sunpy sample data (AIA and RHESSI). Everything worked as expected.

Checked the metadata of the sunpy sample RHESSI data and compared with the Hinode XRT data I downloaded. The metadata were compared because there was no sunpy sample data for XRT and I wanted to make sure that I had included all the necessary information before attempting any plotting.

I saw that I needed to add a few keywords to the Hinode XRT header, which I did [this is after prepping both AIA and XRT in SSWIDL using "aia_prep" and "xrt_prep", respectively]:

Header keywords that were missing from the Hinode XRT fits and were subsequently included are dsun_obs, hglt_obs, hgln_obs, crln_obs, crlt_obs, rsun_obs, rsun_obs, rsun_ref

Created the respective maps

Selected a set of coordinates for an ROI
xlims_world = [797, 905]*u.arcsec
ylims_world = [-468, -356]*u.arcsec

Bottom Left and Top Right coords for each of the maps
blc1 = SkyCoord(xlims_world[0],ylims_world[0], frame=mapBeThick.coordinate_frame)
trc1 = SkyCoord(xlims_world[1],ylims_world[1], frame=mapBeThick.coordinate_frame)
blc2 = SkyCoord(xlims_world[0],ylims_world[0], frame=map193.coordinate_frame)
trc2 = SkyCoord(xlims_world[1],ylims_world[1], frame=map193.coordinate_frame)

SubMaps
subBeThick        = mapBeThick.submap(blc1, top_right=trc1)
sub193                = map193.submap(blc2, top_right=trc2)

Plotted them side-by-side [attached figure "XRT-AIA.png"]

Over-plotted (referring the example from
fig=plt.figure()
ax1 = fig.add_subplot(111, projection=subBeThick)
subBeThick.plot(axes=ax1, cmap='Blues')
sub193.plot(axes=ax1, autoalign=True,alpha=0.5)
plt.show()

The unfortunate-looking plot is attached [ " overplotXRTAIA.png  " ]

The result is the same when the reference frame is changed from subBeThick to sub193
[ "overplotAIAXRT.png" ]

The other issue was these "over-plots" failed to close in the normal way and only responded to a Keyboard Interrupt.

The figures are attached.

I would really, really like some help !
Thanks !

Regards,
Sudha
XRT-AIA.png
overplotAIAXRT.png
overplotXRTAIA.png

Albert Y. Shih

unread,
Mar 22, 2022, 10:17:13 AM3/22/22
to su...@googlegroups.com
Hi, Sudha,
     The reason that the `overplotXRTAIA.png` is "unfortunate-looking" is that you are overplotting imagery at the limb from two different observer locations (SDO and Hinode).  Without applying an additional assumption, the default transformation will discard the AIA data that is off-disk because those pixels are not uniquely defined in the XRT coordinate frame.  The recommended way to retain the off-disk AIA data is to use the `Helioprojective.assume_spherical_screen()` context manager around the plotting code, e.g.:

with Helioprojective.assume_spherical_screen(subBeThick.observer_coordinate, only_off_disk=True):
     <insert plotting code here>

The problem is analogous in `overplotAIAXRT.png` with the off-disk XRT data being discarded because it is not uniquely defined in the AIA coordinate frame.  You may also have an issue in this second plot with your ordering of plotting commands or the setting of alphas.
     I'm also wondering about all of those keywords you set in the XRT header, which are the metadata for the observer location.  (Incidentally, SunPy will ignore `CRLT_OBS`/`CRLN_OBS` if `HGLT_OBS`/`HGLN_OBS` are present.)  If they are missing in the XRT header, how did you know what values to set them to?  Did you copy them over from one of the other spacecraft (e.g., SDO or RHESSI), which of course are not at the same location as Hinode?  I'll note that if you hadn't added those keywords at all, SunPy would have assumed that the observer was at Earth center, which is probably no worse an assumption than using the wrong spacecraft location.
     I don't know why your overplots wouldn't close normally.  My best guess without more information is that the `autoalign` functionality in an interactive plot window can be extremely intensive on the CPU, especially if you don't realize that you are triggering redraws.  You could alternatively try using `reproject_to()` as in some of our other gallery examples, e.g., to create a map that is the XRT data reprojected into the AIA coordinate frame.

Albert

--
You received this message because you are subscribed to the Google Groups "SunPy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sunpy+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sunpy/2dbc48c9-0ae0-4af6-bc9f-306fd167ca97n%40googlegroups.com.

Sudha

unread,
Mar 23, 2022, 5:42:59 AM3/23/22
to SunPy
Hello Albert,

Thanks very much for your explanation and for the recommendation to use "Helioprojective.assume_spherical_screen()".

I considered sub193 as the "base" map or reference and combined "reproject_to()" [ as shown in the example gallery for Reprojecting using a spherical screen ].

with Helioprojective.assume_spherical_screen(sub193.observer_coordinate, only_off_disk=True):
        outBeThick        = subBeThick.reproject_to(sub193.wcs)
        out131               = sub131.reproject_to(sub193.wcs)

Everything comes together beautifully (XRT as contours),
finalOverplot.png
I did not need to use "autoalign" while plotting because of the neat and tidy "reproject_to()"  as you suggested. I haven't yet checked if the option will work without previous issues, now that the maps are projected properly.

On how I set the XRT header keyword values:
from sunpy.sun import constants
from sunpy.coordinates.sun import earth_distance, B0, L0, angular_radius
file1    = 'L1_XRT20110609_035252.0.fits'
fits1    = fits.open(file1)
data1  = fits1[0].data
hdr1    = fits1[0].header
# Writing out the variables for the desired keywords
obs_time    = hdr1['date_obs']
crln_obs     = L0(obs_time).value
crlt_obs      = B0(obs_time).value
rsun_obs    = angular_radius(obs_time).value
rsun_ref      = constants.get('radius').value
dsun_ref     = constants.get('mean distance').value
dsun_obs   = (earth_distance(obs_time)).to_value('m')


Then I set these keywords into the header file and made the map as usual. I'll need to try this procedure with the corresponding RHESSI image and see if it takes.

Thanks !

Regards,
Sudha

Albert Y. Shih

unread,
Mar 23, 2022, 9:13:11 AM3/23/22
to su...@googlegroups.com
Hi, Sudha,
     Great that the overlays are working!
     Regarding your header modifications, you've manually set the observer location to Earth center, which is likely sufficiently accurate for your purposes.  As I noted in my previous email, if you hadn't added those keywords, SunPy would have assumed an observer location of Earth center, so the only difference is that you avoided seeing a SunPy warning saying that it was making that assumption.
     One suggestion: Instead of setting `CRLN_OBS` and `CRLT_OBS`, I recommend you set `HGLN_OBS` (to zero) and `HGLT_OBS` (to `B0(obs_time.value)`).  Not everyone agrees on the calculation of Carrington longitude, but everyone agrees on Stonyhurst longitude.  It makes no difference if you're using only SunPy, but it could make a difference if you take your header values to something other than SunPy.

Albert


Sudha

unread,
Mar 24, 2022, 5:16:01 AM3/24/22
to SunPy
Hi Albert,

Yes, you're right; before I added in the "missing" keywords, the sunpy warning would come up everytime I plotted the XRT map. My main aim then was to figure out the over-plotting issue. Oh yes, sure ! I appreciate your recommendation about the 'hgln_obs' and 'hglt_obs'. I shall study these parameters more and try to understand why and when each one may be used. And explore more about observer locations and subsequent projections.

Thanks so much !

Regards,
Sudha
Reply all
Reply to author
Forward
0 new messages