Guide to Data#
A notebook guide to working with PUNCH FITS files in Python
This Jupyter Notebook presents a comprehensive guide to the analysis and visualization of PUNCH FITS files using Python. The guide leverages a suite of Python libraries, including matplotlib for data visualization, numpy for numerical operations, astropy for FITS file handling and WCS manipulation, sunpy for solar data analysis, and ndcube for multi-dimensional data handling. The notebook serves as a robust resource for researchers and scientists, providing them with the tools
and knowledge to effectively manipulate and interpret PUNCH FITS files. The guide’s step-by-step approach, coupled with code snippets and detailed explanations, ensures a user-friendly experience, promoting the accessibility of PUNCH data for a wider scientific community.
To install these dependencies, use pip install -r requirements.txt in your terminal.
[18]:
# Load libraries
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import numpy as np
from copy import deepcopy
from astropy.io import fits
from astropy.wcs import WCS
import astropy.units as u
from sunpy.map import Map
from ndcube import NDCube
Data loading#
We first specify the file we wish to load. For this notebook, let’s use the included L3 PAM file. Then, we’ll load it two different ways:
With
astropy.ioWith
SunPy Maps
The data are RICE compressed, so we’ll have to look in the second at index 1 (since we’re zero-indexed) HDU for the main data frame. There is also an uncertainty layer in the last HDU that we’ll explore later.
[2]:
filename = 'example_data/PUNCH_L3_PAM_20240620000000.fits'
with fits.open(filename) as hdul:
data = hdul[1].data
header = hdul[1].header
uncertainty = hdul[2].data
First, let’s look at the shape of the data. For this data product, total brightness and polarized brightness are stacked along the first dimension. The uncertainty array corresponds on a pixel-to-pixel basis with the data array.
[3]:
data.shape, uncertainty.shape
[3]:
((2, 4096, 4096), (2, 4096, 4096))
We can look in the header for details about the data. A few things to note:
The header is divided into helpful sections such as “temporal information” and “instrument and spacecraft state.” This aids in navigating the information. These sections are consistent for other products (with the addition of more sections depending on the product).
There are two world coordinate systems (WCS), one for the helio system and one for a celestial right ascension and declination system. They match each other and can be used interchangeably depending on the type of analysis.
The first dimension iterates over total brightness and polarized brightness, in that order. This is indicated by the
OBSLAYR1andOBSLAYR2keywords.
[4]:
header
[4]:
SIMPLE = 'T ' / Conforms to FITS Standard
BITPIX = -32 / Number of bits per pixel
NAXIS = 3 / Number of axes
NAXIS1 = 4096 / Length of the first axis
NAXIS2 = 4096 / Length of the second axis
NAXIS3 = 2 / Length of the third axis
COMMENT ----- FITS Required ----------------------------------------------------
EXTNAME = 'PRIMARY DATA ARRAY' / Name of this binary table extension
LONGSTRN= 'OGIP 1.0' / The OGIP long string convention may be used
COMMENT ----- Documentation, Contact, and Collection Metadata ------------------
DOI = 'https://doi.org/TBD' / Data reference DOI
PROJECT = 'PUNCH '
TITLE = 'PUNCH Level-3 Polarized Low Noise Mosaic'
KEYVOCAB= 'Unified Astronomy Thesaurus Keywords'
KEYWORDS= 'Solar Corona (1483), Solar K Corona (2042), Solar F Corona (1991), &'
CONTINUE 'Solar Coronal Streamers (1486), Solar Coronal Plumes (2039), Solar &'
CONTINUE 'Wind (1534), Fast Solar Wind (1872), Slow Solar Wind (1873), Solar &'
CONTINUE 'Coronal Mass Ejection (310), Heliosphere (711), Polarimetry (1278)'
LICENSE = 'Creative Commons Attribution 4.0 International | CC BY 4.0'
DESCRPTN= 'PUNCH Level-3 data, Composite mosaic in output coordinates'
DOC_URL = 'https://punch.spaceops.swri.org'
COMMENT ----- File Type and Provenance -----------------------------------------
FILENAME= '' / Name of file
LEVEL = '3 ' / Product Level
OBSTYPE = 'Polarized low noise mosaic' / Plain text observation
TYPECODE= 'PA ' / Observation product type code
OBSCODE = 'M ' / Observatory spacecraft code
PIPEVRSN= '' / PUNCHPipe software version number
FILE_RAW= '' / Raw telemetry filename
ORIGIN = 'SwRI ' / Institution responsible for creating the file
COMMENT ----- Temporal Information ---------------------------------------------
TIMESYS = 'UTC ' / Principal time system
DATE-BEG= '2024-06-20T00:00:00.000' / UTC time observation
DATE-OBS= '2024-06-20T00:00:00.000' / UTC reference time
DATE-AVG= '2024-06-20T00:16:00.000' / UTC reference time
DATE-END= '2024-06-20T00:32:00.000' / UTC time of observation end
DATE = '2024-06-20T12:32:00.000' / UTC file generation date and time
COMMENT ----- Instrument and Spacecraft State ----------------------------------
WAVELNTH= 530 / [nm] average peak response
WAVEUNIT= 'nanometer' / Unit of observation measurement
OBS-MODE= 'Polar_BpB' / Image Mode (Unpolarized, Polar_MZP, Polar_BpB)
OBSLAYR1= 'Polar_B ' / Image Mode for first datacube layer
OBSLAYR2= 'Polar_pB' / Image Mode for second datacube layer
INSTRUME= 'WFI+NFI Mosaic' / Instrument name
TELESCOP= 'PUNCH 1-2-3-4' / Satellite name
OBSRVTRY= 'PUNCH ' / Observatory name
OBJECT = 'Heliosphere white light' / Object observed
COMMENT ----- World Coordinate System ------------------------------------------
WCSAXES = 3 / Number of coordinate axes
CRPIX1 = 2047.5 / Pixel coordinate of reference point
CRPIX2 = 2047.5 / Pixel coordinate of reference point
CRPIX3 = 0.0 / Pixel coordinate of reference point
PC1_1 = 1.0 / Coordinate transformation matrix element
PC1_2 = 0.0 / Coordinate transformation matrix element
PC1_3 = 0.0 / Coordinate transformation matrix element
PC2_1 = 0.0 / Coordinate transformation matrix element
PC2_2 = 1.0 / Coordinate transformation matrix element
PC2_3 = 0.0 / Coordinate transformation matrix element
PC3_1 = 0.0 / Coordinate transformation matrix element
PC3_2 = 0.0 / Coordinate transformation matrix element
PC3_3 = 1.0 / Coordinate transformation matrix element
CDELT1 = 0.0225 / [deg] Coordinate increment at reference point
CDELT2 = 0.0225 / [deg] Coordinate increment at reference point
CDELT3 = 1.0 / Coordinate increment at reference point
CUNIT1 = 'deg ' / Units of coordinate increment and value
CUNIT2 = 'deg ' / Units of coordinate increment and value
CUNIT3 = '' / Units of coordinate increment and value
CTYPE1 = 'HPLN-ARC' / Coordinate type codezenithal/azimuthal equidist
CTYPE2 = 'HPLT-ARC' / Coordinate type codezenithal/azimuthal equidist
CTYPE3 = 'STOKES ' / Coordinate type
CRVAL1 = 0.0 / [deg] Coordinate value at reference point
CRVAL2 = 0.0 / [deg] Coordinate value at reference point
CRVAL3 = 0.0 / Coordinate value at reference point
LONPOLE = 180.0 / [deg] Native longitude of celestial pole
LATPOLE = 0.0 / [deg] Native latitude of celestial pole
CNAME3 = 'STOKES ' / Axis name for labelling purposes
MJDREF = 0.0 / [d] MJD of fiducial time
WCSAXESA= 3 / Number of coordinate axes
CRPIX1A = 2047.5 / Pixel coordinate of reference point
CRPIX2A = 2047.5 / Pixel coordinate of reference point
CRPIX3A = 0.0 / Pixel coordinate of reference point
PC1_1A = 0.9912781666299199 / Coordinate transformation matrix element
PC1_2A = 0.13178617667579878 / Coordinate transformation matrix element
PC1_3A = 0.0 / Coordinate transformation matrix element
PC2_1A = -0.13178617667579878 / Coordinate transformation matrix element
PC2_2A = 0.9912781666299199 / Coordinate transformation matrix element
PC2_3A = 0.0 / Coordinate transformation matrix element
PC3_1A = 0.0 / Coordinate transformation matrix element
PC3_2A = 0.0 / Coordinate transformation matrix element
PC3_3A = 1.0 / Coordinate transformation matrix element
CDELT1A = -0.0225 / [deg] Coordinate increment at reference point
CDELT2A = 0.0225 / [deg] Coordinate increment at reference point
CDELT3A = 1.0 / Coordinate increment at reference point
CUNIT1A = 'deg ' / Units of coordinate increment and value
CUNIT2A = 'deg ' / Units of coordinate increment and value
CUNIT3A = ' ' / Units of coordinate increment and value
CTYPE1A = 'RA---ARC' / Coordinate type codezenithal/azimuthal equidist
CTYPE2A = 'DEC--ARC' / Coordinate type codezenithal/azimuthal equidist
CTYPE3A = 'STOKES ' / Coordinate type
CRVAL1A = 88.72529939915754 / [deg] Coordinate value at reference point
CRVAL2A = 23.43076999917664 / [deg] Coordinate value at reference point
CRVAL3A = 0.0 / Coordinate value at reference point
LONPOLEA= 180.0 / [deg] Native longitude of celestial pole
LATPOLEA= 0.0 / [deg] Native latitude of celestial pole
CNAME3A = 'STOKES ' / Axis name for labelling purposes
MJDREFA = 0.0 / [d] MJD of fiducial time
COMMENT ----- Image Statistics and Properties ----------------------------------
BUNIT = 'Mean Solar Brightness' / Units of observation
BSUN_DEF= 20090000.0 / [W/m2/sr] Mean Solar Brightness
DATAZER = 9214166 / number of pixels=0
DATASAT = 0 / number of saturated pixels
DSATVAL = 65535.0 / saturation value in calibrated units
DATAAVG = 4.73192923104541E-14 / average non-zero pixel value
DATAMDN = 2.73807889474573E-15 / median non-zero pixel value
DATASIG = 5.05454673341981E-13 / standard deviation of non-zero pixels
DATAP01 = 6.96018566522324E-16 / 1% of non-zero pixels are LE this value
DATAP10 = 1.07249935352800E-15 / 10% of non-zero pixels are LE this value
DATAP25 = 1.53683609188879E-15 / 25% of non-zero pixels are LE this value
DATAP50 = 2.73807889474573E-15 / 50% of non-zero pixels are LE this value
DATAP75 = 7.08962384038760E-15 / 75% of non-zero pixels are LE this value
DATAP90 = 2.67974000872335E-14 / 90% of non-zero pixels are LE this value
DATAP95 = 7.37835153760170E-14 / 95% of non-zero pixels are LE this value
DATAP98 = 2.76662996982410E-13 / 98% of non-zero pixels are LE this value
DATAP99 = 7.17754143880061E-13 / 99% of non-zero pixels are LE this value
DATAMIN = 0.0 / minimum valid physical value
DATAMAX = 3.78083328533840E-11 / maximum valid physical value
COMMENT ----- Solar Reference Data ---------------------------------------------
RSUN_ARC= / [arcsec] photospheric solar radius
RSUN_REF= 695700000.0 / [m] assumed physical solar radius
SOLAR_EP= -7.424748309496262 / [deg] S/C ecliptic North to solar North angle
CAR_ROT = 2285.6359047394517 / Carrington Rotation number
COMMENT ----- Spacecraft Location & Environment --------------------------------
GEOD_LAT= 0.0 / [deg] S/C sub-point geodetic latitude
GEOD_LON= 0.0 / [deg] S/C sub-point longitude
GEOD_ALT= 0.0 / [m] S/C WGS84 altitude
HGLT_OBS= 1.6391084786225854 / [deg] S/C heliographic latitude
HGLN_OBS= 0.0 / [deg] S/C heliographic longitude (B0)
CRLT_OBS= 1.6391084786225854 / [deg] S/C Carrington latitude (B0)
CRLN_OBS= 131.07429379735413 / [deg] S/C Carrington longitude
DSUN_OBS= 152011862324.1987 / [m] S/C distance from Sun
HEEX_OBS= 152011862324.19867 / [m] S/C Heliocentric Earth Ecliptic X
HEEY_OBS= 2.20584297521290E-05 / [m] S/C Heliocentric Earth Ecliptic Y
HEEZ_OBS= 0.2938932294216845 / [m] S/C Heliocentric Earth Ecliptic Z
HCIX_OBS= -148073749483.16226 / [m] S/C Heliocentric Inertial X
HCIY_OBS= -34100802021.633644 / [m] S/C Heliocentric Inertial Y
HCIZ_OBS= 4348137848.573463 / [m] S/C Heliocentric Inertial Z
HEQX_OBS= 151949662666.6902 / [m] S/C Heliocentric Earth Equatorial X
HEQY_OBS= 0.0 / [m] S/C Heliocentric Earth Equatorial Y
HEQZ_OBS= 4348137848.5734625 / [m] S/C Heliocentric Earth Equatorial Z
COMMENT ----- Fixity -----------------------------------------------------------
CHECKSUM= '1cGK2cGI1cGI1cGI' / HDU checksum updated 2024-06-14T14:59:10
DATASUM = '3773111760' / data unit checksum updated 2024-06-14T14:59:10
COMMENT ----- History ----------------------------------------------------------
HISTORY Records of processing from pipeline
We can create Astropy WCS objects from the header. Since the data are three dimensional, the WCS will be three dimensions. Luckily, we can slice WCSes.
[5]:
data_wcs = WCS(header)
data_wcs
[5]:
WCS Keywords
Number of WCS axes: 3
CTYPE : 'HPLN-ARC' 'HPLT-ARC' 'STOKES'
CRVAL : 0.0 0.0 0.0
CRPIX : 2047.5 2047.5 0.0
PC1_1 PC1_2 PC1_3 : 1.0 0.0 0.0
PC2_1 PC2_2 PC2_3 : 0.0 1.0 0.0
PC3_1 PC3_2 PC3_3 : 0.0 0.0 1.0
CDELT : 0.0225 0.0225 1.0
NAXIS : 4096 4096 2
We can also access the celestial WCS by setting key='A'. “
[6]:
data_wcs_celestial = WCS(header, key='A')
data_wcs_celestial
[6]:
WCS Keywords
Number of WCS axes: 3
CTYPE : 'RA---ARC' 'DEC--ARC' 'STOKES'
CRVAL : 88.72529939915754 23.43076999917664 0.0
CRPIX : 2047.5 2047.5 0.0
PC1_1 PC1_2 PC1_3 : 0.9912781666299199 0.13178617667579878 0.0
PC2_1 PC2_2 PC2_3 : -0.13178617667579878 0.9912781666299199 0.0
PC3_1 PC3_2 PC3_3 : 0.0 0.0 1.0
CDELT : -0.0225 0.0225 1.0
NAXIS : 4096 4096 2
We can also use SunPy Maps to load the data. However, since Maps only support 2D images, it will warn us that it’s truncating a dimension of the data and only showing the total brightness. We get a list of Maps where the first entry is the data and the second is the uncertainty.
[7]:
data_map = Map(filename)
data_map[0]
WARNING: SunpyUserWarning: This file contains more than 2 dimensions. Data will be truncated to the first two dimensions. [sunpy.map.mapbase]
WARNING: SunpyMetadataWarning: Could not parse unit string "Mean Solar Brightness" as a valid FITS unit.
See https://docs.sunpy.org/en/stable/how_to/fix_map_metadata.html for how to fix metadata before loading it with sunpy.map.Map.
See https://fits.gsfc.nasa.gov/fits_standard.html for the FITS unit standards. [sunpy.map.mapbase]
[7]:
<sunpy.map.mapbase.GenericMap object at 0x12520dc70>
|
Image colormap uses histogram equalization
Click on the image to toggle between units 
|
||||||||||||||||||||||||
|
If you wish to make a map of the polarized brightness, we can also do that manually.
[8]:
data_map = Map(data[1], header)
data_map
WARNING: SunpyMetadataWarning: Could not parse unit string "Mean Solar Brightness" as a valid FITS unit.
See https://docs.sunpy.org/en/stable/how_to/fix_map_metadata.html for how to fix metadata before loading it with sunpy.map.Map.
See https://fits.gsfc.nasa.gov/fits_standard.html for the FITS unit standards. [sunpy.map.mapbase]
[8]:
<sunpy.map.mapbase.GenericMap object at 0x120b75a30>
|
Image colormap uses histogram equalization
Click on the image to toggle between units 
|
||||||||||||||||||||||||
|
Plotting data#
There are many ways to plot the data. We are working on more visualization tools that will be available with the release of the punchbowl package. First, we can plot the SunPy map.
[9]:
data_map.plot(norm='log')
[9]:
<matplotlib.image.AxesImage at 0x30ffefb00>
We can also use Matplotlib directly to modify the plot however we want. Notice that we set the projection so it is coordinate aware. Remember we have to slice the WCS as data_wcs[0] instead of data_wcs because our data are 3D but we want to view only the spatial dimensions.
[10]:
# Plot the data using matplotlib manually
fig, ax = plt.subplots(figsize=(9.5, 7.5), subplot_kw={'projection':data_wcs[0]})
im = ax.imshow(data[1], cmap='Greys_r', norm=LogNorm(vmin=1.77e-15, vmax=3.7e-11))
# set up the axis labels
lon, lat = ax.coords
lat.set_ticks(np.arange(-90, 90, 15) * u.degree)
lon.set_ticks(np.arange(-180, 180, 15) * u.degree)
lat.set_major_formatter('dd')
lon.set_major_formatter('dd')
ax.set_facecolor('black')
ax.coords.grid(color='white', alpha=.25, ls='dotted')
ax.set_xlabel("Helioprojective longitude")
ax.set_ylabel("Helioprojective latitude")
ax.set_title('Synthetic PUNCH Polarized Brightness')
fig.colorbar(im, ax=ax, label='Mean Solar Brightness')
[10]:
<matplotlib.colorbar.Colorbar at 0x310af9910>
Since we set the projection, we could also use the celestial WCS instead.
[11]:
# Plot the data using matplotlib manually, using the celestial frame
fig, ax = plt.subplots(figsize=(9.5, 7.5), subplot_kw={'projection':data_wcs_celestial[0]})
im = ax.imshow(data[0], cmap='Greys_r', norm=LogNorm(vmin=1.77e-15, vmax=3.7e-11))
# set up the axis labels
lon, lat = ax.coords
lat.set_ticks(np.arange(-90, 90, 15) * u.degree)
lon.set_ticks(np.arange(-180, 180, 15) * u.degree)
lat.set_major_formatter('dd')
lon.set_major_formatter('dd')
ax.set_facecolor('black')
ax.coords.grid(color='white', alpha=.25, ls='dotted')
ax.set_xlabel("RA")
ax.set_ylabel("DEC")
ax.set_title('Synthetic PUNCH Polarized Brightness')
fig.colorbar(im, ax=ax, label='Mean Solar Brightness')
[11]:
<matplotlib.colorbar.Colorbar at 0x310bdbce0>
Reprojecting the data#
What if we want want to manipulate the images? We can use built-in tools of SunPy Maps. For example, we can reproject.
[12]:
# What if we wanted to reproject the data to a new arbitrary coordinate frame?
# First, define a new target WCS
new_wcs = deepcopy(data_map.wcs)
new_wcs.wcs.ctype = 'HPLN-ARC', 'HPLT-ARC'
new_wcs.wcs.cunit = 'deg', 'deg'
new_wcs.array_shape = 1024, 1024
new_wcs.wcs.crpix = 512, 512
new_wcs.wcs.crval = 10, 10
new_wcs.wcs.cdelt = 0.0225, 0.0225
new_wcs.wcs.pc = (0.66,-0.33), (0.33,0.66)
# We can then reproject using SunPy maps
new_map = data_map.reproject_to(new_wcs)
new_map
[12]:
<sunpy.map.mapbase.GenericMap object at 0x310caaae0>
|
Image colormap uses histogram equalization
Click on the image to toggle between units 
|
||||||||||||||||||||||||
|
Using NDCube#
So far, we’ve only talked about directly manipulating the data or using SunPy Maps. However, we encourage you to utilize NDCubes. We use them in the data reduction pipeline too. NDCube is a SunPy package for generalized n-dimensional data. We can use our manually loaded data to make an NDCube with all of the data in one place. There’s even a way to bundle uncertainty directly into the cube.
[13]:
data_ndcube = NDCube(data, wcs=data_wcs, meta=header, uncertainty=uncertainty)
data_ndcube
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]
[13]:
<ndcube.ndcube.NDCube object at 0x310cabb60>
NDCube
------
Dimensions: [2.000e+00 4.096e+03 4.096e+03] pix
Physical Types of Axes: [('phys.polarization.stokes',), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat'), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat')]
Unit: None
Data Type: float32
We can then use all the great NDCube functionality! Much like SunPy Maps we can reproject even. Note that as of writing this guide, the uncertainty is dropped when NDCube reprojects data.
[14]:
# make a target WCS
new_wcs = deepcopy(data_wcs)
new_wcs.wcs.ctype = 'HPLN-ARC', 'HPLT-ARC', 'STOKES'
new_wcs.wcs.cunit = 'deg', 'deg', ''
new_wcs.array_shape = 2, 1024, 1024
new_wcs.wcs.crpix = 0, 0, 0
new_wcs.wcs.crval = 0, 0, 0
new_wcs.wcs.cdelt = 0.0225, 0.0225, 1.0
new_cube = data_ndcube.reproject_to(new_wcs)
new_cube
[14]:
<ndcube.ndcube.NDCube object at 0x3106c8260>
NDCube
------
Dimensions: [ 2. 1024. 1024.] pix
Physical Types of Axes: [('phys.polarization.stokes',), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat'), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat')]
Unit: None
Data Type: float64
NDCubes have a plotter on them too. Notice we even get an interactive slider to manipulate the polarization axis.
[15]:
new_cube.plot(interpolation='None', norm=LogNorm(vmin=1.77E-15, vmax=3.7E-11), cmap='Greys_r')
[15]:
<mpl_animators.wcs.ArrayAnimatorWCS at 0x310ae6b40>
We encourage you to explore NDCube more. There are other helpful things like cropping that you can do.
Uncertainty#
Level 1, 2, and 3 PUNCH products come with an uncertainty estimate. This uncertainty is expressed as the fractional uncertainty ranging from 0 to 1. To get the absolute uncertainty, simply multiply the fractional uncertainty by the data value.
There’s an extra step though. Uncertainty is encoded as an integer in the FITS file to save space. When we inspect the uncertainty, we see it ranges from 0 to 255. Here, 255 corresponds to a fractional uncertainty of 1.0 or complete uncertainty.
[16]:
print(uncertainty.min(), uncertainty.max()) # values range from 0 to 255
absolute_uncertainty = uncertainty.astype(float) / 255 * data
print(absolute_uncertainty.min(), absolute_uncertainty.max()) # now the uncertainty is in the same units as the data
1 255
0.0 3.780833285338403e-11
Let’s plot the data and uncertainty to get a feel for it. As mentioned earlier, analysis and visualization tools that make this easier will be released with the punchbowl package.
[17]:
vmin, vmax = 1.77E-15, 3.7E-11
fig, axs = plt.subplots(figsize=(14, 6), ncols=2, subplot_kw={'projection':data_wcs[0]})
im = axs[0].imshow(data[0], norm=LogNorm(vmin=vmin, vmax=vmax), interpolation='None', cmap='Greys_r')
axs[0].set_title("Data")
axs[1].imshow(absolute_uncertainty[0], norm=LogNorm(vmin=vmin, vmax=vmax), interpolation='None', cmap='Greys_r')
axs[1].set_title("Uncertainty")
fig.colorbar(im, ax=axs, label='Mean Solar Brightness')
# set up the axis labels
for ax in axs:
lon, lat = ax.coords
lat.set_ticks(np.arange(-90, 90, 15) * u.degree)
lon.set_ticks(np.arange(-180, 180, 15) * u.degree)
lat.set_major_formatter('dd')
lon.set_major_formatter('dd')
ax.set_facecolor('black')
ax.coords.grid(color='white', alpha=.25, ls='dotted')
ax.set_xlabel("Helioprojective longitude")
ax.set_ylabel("Helioprojective latitude")
