This example shows how you can work with FITS data easily and generate plots from the raw data.
Download a sample FITS file from the Sloan Digital Sky Survey (SDSS) data server and save it to a local file:
In[]:=
file=URLDownload["https://data.sdss.org/sas/dr16/sdss/spectro/redux/26/spectra/1323/spec-1323-52797-0012.fits","~/example.fits"];
Examine the general information for this dataset. Note how the date, right ascension, and declination are included:
In[]:=
information=Dataset[Import[file,{"FITS","MetaInformation"}]][1,"GENERAL INFORMATION"]
Out[]=
SIMPLE
True
conforms to FITS standard
BITPIX
8
array data type
NAXIS
0
number of array dimensions
EXTEND
True
—
TAI
4561539252.
1st row - Number of seconds since Nov 17 1858
RA
208.774
1st row - Right ascension of telescope boresigh
DEC
55.1213
1st row - Declination of telescope boresight (d
EQUINOX
2000.0
—
RADECSYS
FK5
—
TAIHMS
07:24:57.65
1st row - TAI time (HH:MM:SS.SS) (TAI-UT = appr
TIMESYS
tai
TAI, not UTC
MJD
52797
MJD of observation
MJDLIST
52792 52797
—
VERSION
NOCVS:v3_113_1:/p/iop/v3_113_1
version of IOP
CAMVER
SPEC1 v4_8
Camera code version
OBSERVER
ank
—
OBSCOMM
science
—
TELESCOP
SDSS 2.5-M
Sloan Digital Sky Survey
BESTEXP
21346
—
NEXP_B1
3
b1 camera number of exposures
rows 1–20 of 139
The FITS data contains multiple spectra. Let’s import one from the second header data unit (HDU):
In[]:=
spectrum=Import[file,{"FITS","Data",2}];
Also get the headers for this data:
In[]:=
headers=Import[file,{"FITS","TableHeaders",2}]
Out[]=
{flux,loglam,ivar,and_mask,or_mask,wdisp,sky,model}
Construct a dataset from this data, using the header names:
In[]:=
ds=Dataset[AssociationThread[headers->Transpose[spectrum]]]
Out[]=
flux
{
…
3827
}
loglam
{
…
3827
}
ivar
{
…
3827
}
and_mask
{
…
3827
}
or_mask
{
…
3827
}
wdisp
{
…
3827
}
sky
{
…
3827
}
model
{
…
3827
}
Convert the “loglam” (the
log
10
wavelengths) column to wavelengths in Angstroms:
In[]:=
wavelengths=QuantityArray[Normal[10^ds["loglam"]],"Angstroms"]
Out[]=
QuantityArray
Dimensions: {3827}
Unit: Angstroms

Compute the measured flux:
In[]:=
flux=QuantityArray[Normal[10^-17*ds["flux"]],"Ergs"/("Centimeters"^2"Seconds""Angstroms")]
Out[]=
QuantityArray
Dimensions: {3827}
Unit:
Ergs
Angstroms
2
Centimeters
Seconds

The conversion factor comes from the metadata information:
In[]:=
information["BUNIT",1]
Out[]=
1E-17 erg/cm^2/s/Ang
Plot the wavelengths versus the flux:
In[]:=
ListLinePlot[Transpose[{wavelengths,flux}],​​PlotRange->All,Frame->True,​​FrameLabel->{QuantityForm["Angstroms","Abbreviation"],QuantityForm["Ergs"/("Centimeters"^2"Seconds""Angstroms"),"Abbreviation"]}]
Out[]=

CITE THIS NOTEBOOK

Plotting an Emission Line Spectrum from SDSS data​
by Arnoud Buzing​
Wolfram Community, STAFF PICKS, January 23, 2024
​https://community.wolfram.com/groups/-/m/t/3107953