Note
Go to the end to download the full example code. or to run this example in your browser via Binder
Guide to QuickPUNCH Data#
A notebook guide to working with QuickPUNCH data in Python
Load libraries
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from ndcube import NDCube
from sunpy.map import Map
from punchbowl.data.sample import QUICKPUNCH_NQN, QUICKPUNCH_WQM
Files Downloaded: 0%| | 0/1 [00:00<?, ?file/s]
PUNCH_LQ_NQN_20230704022400.fits: 0%| | 0.00/1.77M [00:00<?, ?B/s]
PUNCH_LQ_NQN_20230704022400.fits: 6%|▌ | 105k/1.77M [00:00<00:01, 999kB/s]
PUNCH_LQ_NQN_20230704022400.fits: 61%|██████ | 1.08M/1.77M [00:00<00:00, 5.89MB/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00, 2.83file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00, 2.59file/s]
Files Downloaded: 0%| | 0/1 [00:00<?, ?file/s]
PUNCH_LQ_WQM_20230704012000.fits: 0%| | 0.00/1.76M [00:00<?, ?B/s]
PUNCH_LQ_WQM_20230704012000.fits: 5%|▌ | 96.9k/1.76M [00:00<00:01, 965kB/s]
PUNCH_LQ_WQM_20230704012000.fits: 60%|█████▉ | 1.06M/1.76M [00:00<00:00, 5.66MB/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00, 2.89file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00, 2.64file/s]
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(QUICKPUNCH_WQM) 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(QUICKPUNCH_NQN) 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: /home/docs/checkouts/readthedocs.org/user_builds/punchbowl/checkouts/0.0.14/examples/PUNCH_SAMPLE/PUNCH_LQ_WQM_20230704012000.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 4 ()
1 COMPRESSED_IMAGE 1 CompImageHDU 153 (1024, 1024) float32
2 COMPRESSED_IMAGE 1 CompImageHDU 9 (1024, 1024) uint8 (rescales to int8)
NFI QuickPUNCH HDU List:
Filename: /home/docs/checkouts/readthedocs.org/user_builds/punchbowl/checkouts/0.0.14/examples/PUNCH_SAMPLE/PUNCH_LQ_NQN_20230704022400.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 4 ()
1 COMPRESSED_IMAGE 1 CompImageHDU 154 (1024, 1024) float32
2 COMPRESSED_IMAGE 1 CompImageHDU 9 (1024, 1024) uint8 (rescales to int8)
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)
The corresponding headers can be queried as AstroPy header objects
wfi_qp_header["DATE-OBS"]
'2023-07-04T01:20:00.00'
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)
/home/docs/checkouts/readthedocs.org/user_builds/punchbowl/envs/0.0.14/lib/python3.12/site-packages/astropy/wcs/wcs.py:537: FITSFixedWarning: CROTA = 0.0 / Rotation with respect to reference angle
keyword looks very much like CROTAn but isn't.
wcsprm = _wcs.Wcsprm(
/home/docs/checkouts/readthedocs.org/user_builds/punchbowl/envs/0.0.14/lib/python3.12/site-packages/astropy/wcs/wcs.py:537: FITSFixedWarning: CROTA_A = 0.0 / Rotation with respect to reference angle
keyword looks very much like CROTAn but isn't.
wcsprm = _wcs.Wcsprm(
/home/docs/checkouts/readthedocs.org/user_builds/punchbowl/envs/0.0.14/lib/python3.12/site-packages/astropy/wcs/wcs.py:805: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 60129.055556 from DATE-OBS.
Set MJD-BEG to 60129.055556 from DATE-BEG.
Set MJD-AVG to 60129.058333 from DATE-AVG.
Set MJD-END to 60129.061111 from DATE-END'.
warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/punchbowl/envs/0.0.14/lib/python3.12/site-packages/astropy/wcs/wcs.py:805: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 60129.100000 from DATE-OBS.
Set MJD-BEG to 60129.100000 from DATE-BEG.
Set MJD-AVG to 60129.102778 from DATE-AVG.
Set MJD-END to 60129.105556 from DATE-END'.
warnings.warn(
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)
Display this SunPy Map object
wfi_qp_data_map
<sunpy.map.mapbase.GenericMap object at 0x7232538deff0>
|
Image colormap uses histogram equalization
Click on the image to toggle between units 
|
||||||||||||||||||||||||
|
Display this SunPy Map object
nfi_qp_data_map
<sunpy.map.mapbase.GenericMap object at 0x723276d46900>
|
Image colormap uses histogram equalization
Click on the image to toggle between units 
|
||||||||||||||||||||||||
|
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)
Take a quick look at these NDCube objects
wfi_qp_data_ndcube, nfi_qp_data_ndcube
(<ndcube.ndcube.NDCube object at 0x723250c34e90>
NDCube
------
Shape: (1024, 1024)
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 0x723251017200>
NDCube
------
Shape: (1024, 1024)
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)
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")
/home/docs/checkouts/readthedocs.org/user_builds/punchbowl/checkouts/0.0.14/examples/data_guide_quickpunch.py:92: RuntimeWarning: divide by zero encountered in log
plt.imshow(np.log(wfi_qp_data), cmap="Greys_r", vmin=-16, vmax=0)
Text(0.5, 1.0, 'QuickPUNCH Mosaic total brightness - 2023-07-04T01:20:00.00UT')
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")
/home/docs/checkouts/readthedocs.org/user_builds/punchbowl/checkouts/0.0.14/examples/data_guide_quickpunch.py:110: RuntimeWarning: divide by zero encountered in log
plt.imshow(np.log(nfi_qp_data), cmap="Greys_r", vmin=-16, vmax=0)
/home/docs/checkouts/readthedocs.org/user_builds/punchbowl/checkouts/0.0.14/examples/data_guide_quickpunch.py:110: RuntimeWarning: invalid value encountered in log
plt.imshow(np.log(nfi_qp_data), cmap="Greys_r", vmin=-16, vmax=0)
Text(0.5, 1.0, 'QuickPUNCH NFI total brightness - 2023-07-04T02:24:00.00UT')
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 retrieved using astropy.io.fits, if needed, using the disable_image_compression keyword.
with fits.open(QUICKPUNCH_WQM, disable_image_compression=True) as hdul:
header_compression = hdul[1].header
header_compression
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 = 1682295 / 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(1999)' / data format of field: variable length array
TTYPE2 = 'GZIP_COMPRESSED_DATA' / label for field 2
TFORM2 = '1PB(3852)' / 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_20230704012000.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-04T01:20:00.00' / UTC time of observation start
DATE-OBS= '2023-07-04T01:20:00.00' / UTC reference time
DATE-AVG= '2023-07-04T01:24:00.00' / Mean UTC observation time
DATE-END= '2023-07-04T01:28:00.00' / UTC time of observation end
DATE = '2023-07-04T13:28: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-ARC' / Coordinate type codezenithal/azimuthal equidist
CTYPE2 = 'HPLT-ARC' / Coordinate type codezenithal/azimuthal equidist
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 = 2048.0 / Pixel coordinate of reference point
CRPIX2 = 2048.0 / 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.0225 / [deg] Coordinate increment at reference point
CDELT2 = 0.0225 / [deg] Coordinate increment at reference point
PC1_1 = 0.9944226355 / WCS coordinate description matrix
PC1_2 = -0.0572686183 / WCS coordinate description matrix
PC2_1 = 0.05721285161 / WCS coordinate description matrix
PC2_2 = 0.994136111 / 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= 152091235598.33926 / [m] Obs distance from Sun
CRLT_OBS= 0.0 / [deg] Obs Carrington latitude
CRLN_OBS= 95.8156043715748 / [deg] Obs Carrington longitude
HEEX_OBS= 694649626.5115166 / [m] Obs Heliocentric Earth Ecliptic X
HEEY_OBS= -4350099.9957544515 / [m] Obs Heliocentric Earth Ecliptic Y
HEEZ_OBS= -37966603.97556185 / [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
Total running time of the script: (0 minutes 2.660 seconds)
