QuickPUNCH data guide#
A notebook guide to working with QuickPUNCH data in Python
[1]:
# Load libraries
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from astropy.io.fits import header
import astropy.units as u
from sunpy.map import Map
from ndcube import NDCube
[2]:
# Specify data filepath
wfi_qp_filename = './data/PUNCH_L2_WQM_20230704000000.fits'
nfi_qp_filename = './data/PUNCH_L2_NQN_20230704000000.fits'
[3]:
# Open the HDU list, and read out the appropriate data
# As the data is RICE compressed, the *second* HDU contains the main data frame
# The third HDU contains a corresponding uncertainty array
with fits.open(wfi_qp_filename) as hdul:
print('WFI QuickPUNCH HDU List:')
hdul.info()
wfi_qp_data = hdul[1].data
wfi_qp_header = hdul[1].header
wfi_qp_uncertainty = hdul[2].data
with fits.open(nfi_qp_filename) as hdul:
print('NFI QuickPUNCH HDU List:')
hdul.info()
nfi_qp_data = hdul[1].data
nfi_qp_header = hdul[1].header
nfi_qp_uncertainty = hdul[2].data
WFI QuickPUNCH HDU List:
Filename: ./data/PUNCH_L2_WQM_20230704000000.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 4 ()
1 COMPRESSED_IMAGE 1 CompImageHDU 156 (1024, 1024) float32
2 COMPRESSED_IMAGE 1 CompImageHDU 9 (1024, 1024) uint8
NFI QuickPUNCH HDU List:
Filename: ./data/PUNCH_L2_NQN_20230704000000.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 4 ()
1 COMPRESSED_IMAGE 1 CompImageHDU 157 (1024, 1024) float32
2 COMPRESSED_IMAGE 1 CompImageHDU 9 (1024, 1024) uint8
[4]:
# The primary data arrays are stored as standard ndarrays
# The uncertainty data array has the dimensions as the primary data array
# Both the primary and uncertainty data arrays share the same header, contained in the primary HDU
print('WFI data array size:', wfi_qp_data.shape)
print('WFI uncertainty array size:', wfi_qp_uncertainty.shape)
WFI data array size: (1024, 1024)
WFI uncertainty array size: (1024, 1024)
[5]:
# The corresponding headers can be queried as AstroPy header objects
wfi_qp_header['DATE-OBS']
[5]:
'2023-07-04T00:00:00.00'
[6]:
# The header information can be converted into an AstroPy WCS object
wfi_qp_data_wcs = WCS(wfi_qp_header);
nfi_qp_data_wcs = WCS(nfi_qp_header);
WARNING: FITSFixedWarning: CROTA = 0.0 / Rotation with respect to reference angle
keyword looks very much like CROTAn but isn't. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: CROTA_A = 0.0 / Rotation with respect to reference angle
keyword looks very much like CROTAn but isn't. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 60129.000000 from DATE-OBS.
Set MJD-BEG to 60129.000000 from DATE-BEG.
Set MJD-AVG to 60129.002778 from DATE-AVG.
Set MJD-END to 60129.005556 from DATE-END'. [astropy.wcs.wcs]
[7]:
# Construct a SunPy Map object of out this data
wfi_qp_data_map = Map(wfi_qp_data, wfi_qp_header)
nfi_qp_data_map = Map(nfi_qp_data, nfi_qp_header)
[8]:
# Display this SunPY Map object
wfi_qp_data_map
WARNING: SunpyMetadataWarning: Could not parse unit string "Mean Solar Brightness" as a valid FITS unit.
See https://docs.sunpy.org/en/stable/code_ref/map.html#fixing-map-metadata for how to fix metadata before loading it with sunpy.map.Map.
See https://fits.gsfc.nasa.gov/fits_standard.html forthe FITS unit standards. [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/code_ref/map.html#fixing-map-metadata for how to fix metadata before loading it with sunpy.map.Map.
See https://fits.gsfc.nasa.gov/fits_standard.html forthe FITS unit standards. [sunpy.map.mapbase]
[8]:
<sunpy.map.mapbase.GenericMap object at 0x7fc0097818e0>
|
Image colormap uses histogram equalization
Click on the image to toggle between units 
|
||||||||||||||||||||||||
|
[9]:
# Display this SunPY Map object
nfi_qp_data_map
WARNING: SunpyMetadataWarning: Could not parse unit string "Mean Solar Brightness" as a valid FITS unit.
See https://docs.sunpy.org/en/stable/code_ref/map.html#fixing-map-metadata for how to fix metadata before loading it with sunpy.map.Map.
See https://fits.gsfc.nasa.gov/fits_standard.html forthe FITS unit standards. [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/code_ref/map.html#fixing-map-metadata for how to fix metadata before loading it with sunpy.map.Map.
See https://fits.gsfc.nasa.gov/fits_standard.html forthe FITS unit standards. [sunpy.map.mapbase]
[9]:
<sunpy.map.mapbase.GenericMap object at 0x7fc009781730>
|
Image colormap uses histogram equalization
Click on the image to toggle between units 
|
||||||||||||||||||||||||
|
[10]:
# Construct an NDCube object out of this data
wfi_qp_data_ndcube = NDCube(wfi_qp_data, wcs=wfi_qp_data_wcs)
nfi_qp_data_ndcube = NDCube(nfi_qp_data, wcs=nfi_qp_data_wcs)
[11]:
# Take a quick look at these NDCube objects
wfi_qp_data_ndcube, nfi_qp_data_ndcube
[11]:
(<ndcube.ndcube.NDCube object at 0x7fbfe96fe3a0>
NDCube
------
Dimensions: [1024. 1024.] pix
Physical Types of Axes: [('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat'), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat')]
Unit: None
Data Type: float32,
<ndcube.ndcube.NDCube object at 0x7fbfd92aee50>
NDCube
------
Dimensions: [1024. 1024.] pix
Physical Types of Axes: [('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat'), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat')]
Unit: None
Data Type: float32)
[12]:
# Display this data in a regular plotting environment, using the associated WCS
plt.figure(figsize=(7.5, 7.5))
ax = plt.subplot(111, projection=wfi_qp_data_wcs)
plt.imshow(np.log(wfi_qp_data), cmap='Greys_r', vmin=-16, vmax=0)
lon, lat = ax.coords
lat.set_ticks(np.arange(-90, 90, 5) * u.degree)
lon.set_ticks(np.arange(-180, 180, 5) * u.degree)
lat.set_major_formatter('dd')
lon.set_major_formatter('dd')
ax.set_facecolor('black')
ax.coords.grid(color='white', alpha=.1)
plt.xlabel("Helioprojective longitude")
plt.ylabel("Helioprojective latitude")
plt.scatter(0, 0, s=240, color='k', transform=ax.get_transform('world'));
plt.title('QuickPUNCH Mosaic total brightness - ' + wfi_qp_header['DATE-OBS'] + 'UT')
/var/folders/py/f06qj9p54jzfxlzpr79z7z0wvsr0ty/T/ipykernel_68963/3721481304.py:5: RuntimeWarning: divide by zero encountered in log
plt.imshow(np.log(wfi_qp_data), cmap='Greys_r', vmin=-16, vmax=0)
[12]:
Text(0.5, 1.0, 'QuickPUNCH Mosaic total brightness - 2023-07-04T00:00:00.00UT')
[13]:
# Display this data in a regular plotting environment, using the associated WCS
plt.figure(figsize=(7.5, 7.5))
ax = plt.subplot(111, projection=nfi_qp_data_wcs)
plt.imshow(np.log(nfi_qp_data), cmap='Greys_r', vmin=-16, vmax=0)
lon, lat = ax.coords
lat.set_ticks(np.arange(-90, 90, 5) * u.degree)
lon.set_ticks(np.arange(-180, 180, 5) * u.degree)
lat.set_major_formatter('dd')
lon.set_major_formatter('dd')
ax.set_facecolor('black')
ax.coords.grid(color='white', alpha=.1)
plt.xlabel("Helioprojective longitude")
plt.ylabel("Helioprojective latitude")
plt.scatter(0, 0, s=240, color='k', transform=ax.get_transform('world'));
plt.title('QuickPUNCH NFI total brightness - ' + nfi_qp_header['DATE-OBS'] + 'UT')
/var/folders/py/f06qj9p54jzfxlzpr79z7z0wvsr0ty/T/ipykernel_68963/2980633537.py:5: RuntimeWarning: divide by zero encountered in log
plt.imshow(np.log(nfi_qp_data), cmap='Greys_r', vmin=-16, vmax=0)
[13]:
Text(0.5, 1.0, 'QuickPUNCH NFI total brightness - 2023-07-04T00:00:00.00UT')
[14]:
# Again noting that these files are compressed, additional keywords will be visible when viewing these FITS files outside of Python.
# These keywords relate to the compression implementation, and can be retreived using astropy.io.fits, if needed, using the disable_image_compression keyword.
with fits.open(wfi_qp_filename, disable_image_compression=True) as hdul:
header_compression = hdul[1].header
header_compression
[14]:
XTENSION= 'BINTABLE' / binary table extension
BITPIX = 8 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 32 / width of table in bytes
NAXIS2 = 1024 / number of rows in table
PCOUNT = 1602185 / number of group parameters
GCOUNT = 1 / number of groups
TFIELDS = 4 / number of fields in each row
TTYPE1 = 'COMPRESSED_DATA' / label for field 1
TFORM1 = '1PB(1985)' / data format of field: variable length array
TTYPE2 = 'GZIP_COMPRESSED_DATA' / label for field 2
TFORM2 = '1PB(0) ' / data format of field: variable length array
TTYPE3 = 'ZSCALE ' / label for field 3
TFORM3 = '1D ' / data format of field: 8-byte DOUBLE
TTYPE4 = 'ZZERO ' / label for field 4
TFORM4 = '1D ' / data format of field: 8-byte DOUBLE
ZIMAGE = T / extension contains compressed image
ZSIMPLE = T / Conforms to FITS Standard
ZBITPIX = -32 / Number of bits per pixel
ZNAXIS = 2 / Number of axes
ZNAXIS1 = 1024 / Length of the first axis
ZNAXIS2 = 1024 / Length of the second axis
ZTILE1 = 1024 / size of tiles to be compressed
ZTILE2 = 1 / size of tiles to be compressed
ZCMPTYPE= 'RICE_1 ' / compression algorithm
ZNAME1 = 'BLOCKSIZE' / compression block size
ZVAL1 = 32 / pixels per block
ZNAME2 = 'BYTEPIX ' / bytes per pixel (1, 2, 4, or 8)
ZVAL2 = 4 / bytes per pixel (1, 2, 4, or 8)
ZNAME3 = 'NOISEBIT' / floating point quantization level
ZVAL3 = 16.0 / floating point quantization level
ZQUANTIZ= 'NO_DITHER' / No dithering during quantization
EXTNAME = 'COMPRESSED_IMAGE' / name of this binary table extension
COMMENT ----- FITS Required ----------------------------------------------------
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-2 WFI+NFI Mosaic QuickPUNCH Clear Image'
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'
COMMENT PUNCH Level-2 NFI QuickPUNCH Clear data, calibrated for non-linear
COMMENT photometric correction, deficient pixel removal, data destreaking and
COMMENT cosmic-ray despiking, vignetting and linear flat-field correction,
COMMENT stray-light subtraction, PSF correction, star-field alignment, and
COMMENT data resampling
COMMENT Rice data compression used
COMMENT Bad data flagging and quality marking is included in a secondary
COMMENT uncertainty HDU
COMMENT Documentation is available on the PUNCH website:
COMMENT https://punch.spaceops.swri.org and via the DOI referenced above (TBD)
COMMENT ----- File Type and Provenance -----------------------------------------
FILENAME= 'PUNCH_L2_WQM_20230704000000.fits' / Name of file
LEVEL = 2 / Product Level
OBSTYPE = 'Unpolarized WFI+NFI Mosaic QuickPUNCH image' / Plain text observation
TYPECODE= 'WQ' / Observation product type code
L1PAREN1= 'PUNCH_L1_CL1_20220204134800_v0.fits' / name of parent file
L1PAREN2= 'PUNCH_L1_CL2_20220204134700_v0.fits' / name of parent file
L1PAREN3= 'PUNCH_L1_CL3_20220204134600_v0.fits' / name of parent file
L1PAREN4= 'PUNCH_L1_CL4_20220204135100_v0.fits' / name of parent file
PIPEVRSN= '1.2.3' / PUNCHPipe software version number
ORIGIN = 'SwRI' / Institution responsible for creating the file
COMMENT ----- Temporal Information ---------------------------------------------
TIMESYS = 'UTC' / Principle time system
DATE-BEG= '2023-07-04T00:00:00.00' / UTC time of observation start
DATE-OBS= '2023-07-04T00:00:00.00' / UTC reference time
DATE-AVG= '2023-07-04T00:04:00.00' / Mean UTC observation time
DATE-END= '2023-07-04T00:08:00.00' / UTC time of observation end
DATE = '2023-07-04T12:08:00.00' / UTC file generation date and time
COMMENT ----- Instrument and Spacecraft State ----------------------------------
WAVELNTH= 530 / [nm] Average peak response
WAVEUNIT= 'nanometer' / Unit of observation measurement
FILTER = 'Clear' / Name of filter for observations
OBS-MODE= 'Unpolarized' / Image Mode (Unpolarized, Polar_M, _Z, or _P)
INSTRUME= 'WFI+NFI MOSAIC' / Instrument name
TELESCOP= 'PUNCH01-02-03-04' / Satellite name
OBSRVTRY= 'PUNCH' / Observatory name
OBJECT = 'Heliosphere white light' / Object observed
COMMENT ----- World Coordinate System ------------------------------------------
WCSAXES = 2 / Number of coordinate axes
IMAGEW = 1024 / Image width, in pixels
IMAGEH = 1024 / Image height, in pixels
WCSNAME = 'Helioprojective Zenithal-Azimuthal Perspective' / Coordinate System
CTYPE1 = 'HPLN-AZP' / Coordinate type codezenithal/azimuthal perspect
CTYPE2 = 'HPLT-AZP' / Coordinate type codezenithal/azimuthal perspect
EQUINOX = 2000.0 / Equatorial coordinates definition (yr)
LONPOLE = 180.0 / [deg] Native longitude of celestial pole
LATPOLE = 0.0 / [deg] Native latitude of celestial pole
CRVAL1 = 0.0 / [deg] Coordinate value at reference point
CRVAL2 = 0.0 / [deg] Coordinate value at reference point
CRPIX1 = 511.5 / Pixel coordinate of reference point
CRPIX2 = 511.5 / Pixel coordinate of reference point
CUNIT1 = 'deg' / Units of coordinate increment and value
CUNIT2 = 'deg' / Units of coordinate increment and value
CROTA = 0.0 / Rotation with respect to reference angle
CDELT1 = 0.045 / [deg] Coordinate increment at reference point
CDELT2 = 0.045 / [deg] Coordinate increment at reference point
PC1_1 = 1.0 / WCS coordinate description matrix
PC1_2 = 0.0 / WCS coordinate description matrix
PC2_1 = 0.0 / WCS coordinate description matrix
PC2_2 = 1.0 / WCS coordinate description matrix
WCSNAMEA= 'Right Ascension Declination Zenithal-Azimuthal' / Coordinate System
CTYPE1A = 'RA---AZP' / Right Ascension - Zenithal-Azimuth Perspective
CTYPE2A = 'DEC--AZP' / Declination - Zenithal-Azimuth Perspective
EQUINOXA= 2000.0 / Equatorial coordinates definition (yr)
LONPOLEA= 180.0 / [deg] Native longitude of the celestial pole
LATPOLEA= 0.0 / [deg] Native latitude of the celestial pole
CRVAL1A = 249.3460044 / [arcsec] X-coordinate of reference point
CRVAL2A = -27.25011692 / [arcsec] Y-coordinate of reference point
CRPIX1A = 495.773 / Pixel coordinate of reference point
CRPIX2A = 507.555 / Pixel coordinate of reference point
CUNIT1A = 'deg' / WCS axis X units (theta_x, right ascension)
CUNIT2A = 'deg' / WCS axis Y units (theta_y, declination)
CROTA_A = 0.0 / Rotation with respect to reference angle
CDELT1A = -0.042305 / [arcsec] Average pixel scale along axis 1
CDELT2A = 0.042305 / [arcsec] Average pixel scale along axis 2
PC1_1A = 0.9944226355 / WCS coordinate description matrix
PC1_2A = -0.0572686183 / WCS coordinate description matrix
PC2_1A = 0.05721285161 / WCS coordinate description matrix
PC2_2A = 0.994136111 / WCS coordinate description matrix
PV1_1A = 0.0 / [deg] Native longitude of the reference point
PV1_2A = 90.0 / [deg] Native latitude of the reference point
PV1_3A = 180.0 / [deg] Alias for LONPOLE (has precedence)
PV2_0A = -1.99E-08 / ZPN projection parameters m=0
PV2_1A = 1.005010009 / ZPN projection parameters m=1
PV2_2A = 0.07295829803 / ZPN projection parameters m=2
PV2_3A = 0.2752920091 / ZPN projection parameters m=3
PV2_4A = -0.7018809915 / ZPN projection parameters m=4
PV2_5A = 1.97518003 / ZPN projection parameters m=5
COMMENT ----- Image Statistics and Properties ----------------------------------
BUNIT = 'Mean Solar Brightness' / Units of observation [6.4282E7 W/m2/sr]
DATAZER = 0 / Number of pixels=0
DATASAT = 1 / Number of saturated pixels
DSATVAL = 1.92E-10 / Saturation value in calibrated units
DATAAVG = 3.89E-11 / Average non-zero pixel value
DATAMDN = 2.82E-11 / Median non-zero pixel value
DATASIG = 3.10E-11 / Standard deviation of non-zero pixels
DATAP01 = 9.90E-13 / 1% of non-zero pixels are LE this value
DATAP10 = 1.43E-11 / 10% of non-zero pixels are LE this value
DATAP25 = 1.89E-11 / 25% of non-zero pixels are LE this value
DATAP50 = 2.82E-11 / 50% of non-zero pixels are LE this value
DATAP75 = 4.73E-11 / 75% of non-zero pixels are LE this value
DATAP90 = 8.01E-11 / 90% of non-zero pixels are LE this value
DATAP95 = 1.08E-10 / 95% of non-zero pixels are LE this value
DATAP98 = 1.40E-10 / 98% of non-zero pixels are LE this value
DATAP99 = 1.59E-10 / 99% of non-zero pixels are LE this value
DATAMIN = 23 / Minimum valid physical value
DATAMAX = 1234 / Maximum valid physical value
COMMENT ----- Solar Reference Data ---------------------------------------------
RSUN_ARC= 6276.642451 / [arcsec] Photospheric solar radius
RSUN_REF= 695507968 / [m] Assumed physical solar radius
SOLAR_P0= 0.217279871 / [deg] Angle of geocentric north to solar north
CAR_ROT = 2272 / Carrington Rotation number
COMMENT ----- Spacecraft Location & Environment --------------------------------
GEOD_LAT= 0 / [deg] Obs sub-point geodetic latitude
GEOD_LON= 0 / [deg] Obs sub-point longitude
GEOD_ALT= 0 / [m] Obs WGS84 altitude
HGLT_OBS= 0.0 / [deg] Obs heliographic latitude (B0)
HGLN_OBS= 0.0 / [deg] Obs heliographic longitude
DSUN_OBS= 152091156482.06534 / [m] Obs distance from Sun
CRLT_OBS= 0.0 / [deg] Obs Carrington latitude
CRLN_OBS= 96.55094869704908 / [deg] Obs Carrington longitude
HEEX_OBS= 694653643.4177302 / [m] Obs Heliocentric Earth Ecliptic X
HEEY_OBS= -4343692.222176673 / [m] Obs Heliocentric Earth Ecliptic Y
HEEZ_OBS= -37893772.896766946 / [m] Obs Heliocentric Earth Ecliptic Z
COMMENT ----- FIXITY -----------------------------------------------------------
CHECKSUM= '9hAOGg3O9g9OEg9O' / HDU checksum updated 2022-03-23T12:35:22
DATASUM = '4040810964' / Data unit checksum updated 2022-03-23T12:35:22
ZHECKSUM= '9hAOGg3O9g9OEg9O' / HDU checksum updated 2022-03-23T12:35:22
ZDATASUM= '4040810964' / Data unit checksum updated 2022-03-23T12:35:22
COMMENT ----- HISTORY ----------------------------------------------------------
HISTORY Records of processing from pipeline
[ ]: