# A Beginner’s Guide to Working with Astronomical Data

Markus Pössel
Haus der Astronomie and Max Planck Institute for Astronomy

###### Contents

  1. 1 Introduction
    1. 1.1 Types of data
    2. 1.2 Types of tools
    3. 1.3 Concepts and operations
    4. 1.4 Software/language choices
  2. 2 Data basics: images, spectra, tables
    1. 2.1 Images: Colour, brightness, pixels
    2. 2.2 Images: PSF and noise
    3. 2.3 Images: Noise and flatfielding
    4. 2.4 Images: astronomical information
    5. 2.5 Spectra
    6. 2.6 Data cubes
    7. 2.7 High-level data: catalogues and tables
  3. 3 SAOImage DS9 and astronomical images
    1. 3.1 Loading a Hubble image
    2. 3.2 A first look at the Eagle Nebula M16
    3. 3.3 Coordinates: Navigating the image
    4. 3.4 Meta-Information: the FITS header
    5. 3.5 Making a colour image
    6. 3.6 Catalogs
    7. 3.7 Photometry with regions and statistics
    8. 3.8 Profiles
  4. 4 TOPCAT and table data
    1. 4.1 Opening a table file
    2. 4.2 Making a sky plot
    3. 4.3 Virtual Observatory (VO) services
    4. 4.4 Basic ADQL queries
    5. 4.5 Selections and subsets
    6. 4.6 More on plotting
    7. 4.7 Histograms
    8. 4.8 A quick look at a spectrum
  5. 5 Getting started with Python
    1. 5.1 Installing Python
    2. 5.2 Using Python in Spyder
    3. 5.3 Modules
  6. 6 Basic operations with Python
    1. 6.1 Meet your new versatile calculator
    2. 6.2 Units and constants
    3. 6.3 Random numbers
    4. 6.4 Strings
    5. 6.5 Conditions
    6. 6.6 User-defined functions
    7. 6.7 Timing your code
  7. 7 Taming long data sets: Lists in Python
    1. 7.1 A list of galaxies
    2. 7.2 Doing something element by element
    3. 7.3 Operations involving more than one list
    4. 7.4 Creating lists simultaneously
    5. 7.5 Numpy arrays
    6. 7.6 Variable types, lists, arrays and speed
    7. 7.7 Strings and base n numbers as lists
  8. 8 Basic plotting with Python and Matplotlib
    1. 8.1 Plotting a function
    2. 8.2 Making a plot look better
    3. 8.3 Annotating plots
    4. 8.4 Figure size
    5. 8.5 Scatter plots
    6. 8.6 Fitting data
    7. 8.7 Histograms
    8. 8.8 Saving figures
    9. 8.9 Glueing data sets
  9. 9 Importing table data into Python
    1. 9.1 Opening a FITS table in python
    2. 9.2 Opening an ASCII table in python
    3. 9.3 Accessing astronomical data bases
  10. 10 Astronomical image manipulation with Python
    1. 10.1 FITS files and python
    2. 10.2 Displaying (showing) an image
    3. 10.3 Pixelwise operations
  11. 11 A simple simulation
    1. 11.1 Step-by-step numerical integration: Euler method
    2. 11.2 Numerical errors
    3. 11.3 Velocity Verlet algorithm
    4. 11.4 A simple two-dimensional simulation
  12. 12 Conclusion

##  1 Introduction

Compare professional astronomy today with how it was 50 years ago, and you
will recognise some continuity — but also a number of fundamental changes.
Perhaps the key change is that astronomy has come to rely almost completely on
digital data. Modern telescopes with their CCD cameras produce digital images
and, with the help of suitable dispersive elements, digital astronomical
spectra. An in-depth analysis of a particularly well-studied object will be
able to make use of digital images and spectra taken at different wavelengths
— some taken by ground-based telescopes and some, like extreme ultraviolet or
X rays, which can only be provided by space telescopes.

Furthermore, in-depth studies of selected objects are only part of what modern
astronomy has to offer. We also live in an era of extensive surveys: large-
scale undertakings to photograph, or take spectra of, wider regions of the sky
in a systematic way. These surveys not only produce many images and spectra,
but also extensive catalogues of the objects observed, listing various of
their properties. With such catalogues comes the ability to make statistical
deductions about astronomical objects: If you want to know whether, say,
elliptical galaxies are, in general, brighter than spiral galaxies (to pick an
artificially simple example), you, consult a suitable galaxy catalogue and
look up the brightness values for a large number of elliptical and a large
number of spiral galaxies.

Modern surveys produce considerable amounts of data. For home use, we have
become used to Gigabytes (1 Gigabyte =1000 Megabyte) and Terabytes (TB; 1
Terabyte = 1000 Gigabyte): a DVD holds 4.7 Gigabytes, and hard drives now
routinely hold a Terabyte or more. The ESO Survey Telescope VISTA produces
about 150 TB worth of data per year, and the Large Synoptic Survey Telescope
(LSST) currently under construction is predicted to produce 500 TB worth of
image data per month.

Then, there a increasingly large and detailed simulations. Take the
IllustrisTNG simulation, which follows the evolution of a large portion of the
universe from shortly after the Big Bang to the present. The smallest but most
detailed of the TNG runs, TNG50, follows the fate of a cube that, in the
present universe, has a side-length of 50 Mpc. Within this volume, matter is
represented by 10 billion point particles representing Dark Matter and 10
billion point particles representing gas. The two simulation runs TNG300 and
TNG100 which were made available to the public in December 2018 sum up to more
than one Petabyte of data (PB; 1 Petabyte = 1000 Terabytes).

Increase the amount of data, and at some point it will become impractical to
download a complete data set onto your own computer for analysis. This is
where data base operations become important: the data is stored in dedicated
data centres, and is accessible online; in order to work with the data, you
use the Internet to send specific queries (”Give me the list of all galaxies
on the Southern hemisphere which are brighter than X”). In this way, the only
data you download is the data you specifically need for your research. The
next step is not far off: when even those pre-selected data sets become too
cumbersome to handle, researchers can run their analysis programs remotely on
the dedicated servers where the data is stored. Infrastructure to allow just
this, notably JuPyter notebooks, are becoming increasingly common.

All this implies that digital data analysis skills are part of the essential
skill sets of modern astronomers. Some of the skills needed for a given
research project will be very specific, involving custom software to be used
for a very particular kind of analysis, or custom software to be written by
the researcher herself. These special skills must be learned and honed on the
job. But there are other skills which are more elementary and more general.
Teaching a selection of those skills is the purpose of this text. It was
originally written for interns at Haus der Astronomie in Heidelberg, in
particular for participants of our International Summer Internship
Program111http://www.haus-der-astronomie.de/en/what-we-do/internships/summer-
internship aimed at students in the final years of high school, or for
students who have just finished high school and are about to start college.

The text is meant to give a first introduction to working with astronomical
data. It does not cover the more detailed astronomical use cases, but instead
is meant to help students familiarise themselves with the basic tools needed
for such work, and learn to apply basic techniques and tools that are fairly
universal.

###  1.1 Types of data

When it comes to data from observational astronomy, most data sets we will be
dealing with fall into one of the following categories:

  * •

Image data — in its simplest form, an image is a two-dimensional array of
pixels, where each pixel value denotes a brightness value. In an ordinary
color image, each pixel will have three brightness values, denoting the
contribution from red, green, and blue (RGB). Since astronomers use many
specialist filters beyond these three colors, astronomical ”color” images can
have even more color values per pixel. Astronomical images usually show a
region of the night sky.

  * •

Spectra — simple spectra show us how the energy of the light emitted by an
object is distributed among the different possible wavelengths. Such simple
spectra are one-dimensional: for each wavelength value, we know the
contribution of light from that particular wavelength region.

  * •

Data cubes — think of these as an enhanced version of astronomical images. An
example is a data cube from what is known as integral field spectroscopy
(IFS); such a data cube is like a two-dimensional image, but now each pixel
contains not a brightness value, but a whole spectrum received from the region
of the sky within that pixel. Since we have a one-dimensional spectrum for
each pixel of a two-dimensional image, that gives us in effect a three-
dimensional object: a data cube

  * •

Catalog data — on a higher level of analysis, astronomers make catalogues of
the properties of different types of astronomical objects. A star catalog, for
instance, could list position, proper motion, parallax, brightness (in various
wavelength bands), and effective temperature for each of a specific selection
of stars.

This list is not complete — for instance, in interferometric imaging, when you
are trying to reconstruct an image by combining coherently the measurements of
different telescopes (“aperture synthesis”), your raw data will be time-
stamped data from the single telescopes, and the initial processing will
involve cross-correlating the data between those telescopes. But while the
list is not complete, it should cover the great majority of current
astronomical use cases.

We will take a first look at examples for each data type in section 2. These
different data usually come with meta-data, that is, descriptive information
about the data. Astronomical images typically include information about the
circumstances of when and how the image was taken (what telescope, what time,
what exposure time, what pointing?), and about where in the sky the object in
question is located (in the shape of data allowing the user to relate image
pixels to an astronomical coordinate system).

For simulations, the situation is more diverse, but there are two fundamental
schemes:

  * •

N-body simulations — here, matter is represented by point particles. A point
particle can represent different kinds of objects: it could be a lump of gas,
or a star, or a group of a few $10^{4}$ or $10^{5}$ stars in a galaxy, or a
lump of dark matter. But every particle has a position in space, and as the
simulation runs, the particle positions change.

  * •

Grid-based simulations — here, space is divided into basic cells, for instance
a space-filling set of small cubes. For each cell, basic properties (such as
density, temperature, quantities of the different types of matter present) are
tracked; those values change as the simulation runs. In more complex
simulations, the grid itself can also change in ways that are adapted to
making the simulation more efficient (“adaptive grid”).

###  1.2 Types of tools

When it comes to the tools for working with these various kinds of data, we
can distinguish two types.

  * •

Application software is software written for a specific set of tasks. In
everyday electronic life, an image viewer allows us to inspect image data, an
image manipulation program such as Adobe Photoshop or The Gimp allows us to
change images in specific ways, and Microsoft Office Excel is a common way of
dealing with catalog-data in the form of spread-sheets.

Good application software has the advantage of being comparatively easy to
operate — in line with modern usage, such software offers you a menu structure
for the selection of task, and a graphical and interactive interface. Also,
while application software performs only a limited selection of tasks, good
application software performs those tasks rather well, having been written by
specialists who have a lot of experience with the kind of task in question.

Increasingly, web applications are playing an important role in everyday life
— applications that you run in your web browser, with a grey area between more
complex applications and dynamically generated web pages. Dynamic web content
does play an important role in astronomy, as well: Many researchers in
astronomy will begin their day by looking over the new research articles
posted as e-prints at [https://arxiv.org/list/astro-ph/recent], use NASA’s
Astrophysical Data Service ADS [https://ui.adsabs.harvard.edu/] to look up
specific research articles, find information about specific objects in the
CDS’s SIMBAD astronomical data base [http://simbad.u-strasbg.fr/simbad/] or
the NASA/IPAC extragalactic data base NED [https://ned.ipac.caltech.edu/] —
and, along the rest of us, use Google as a search tool or read up on helpful
answers on sites such as Stack Overflow.

We will use some application software in the following, namely SAOImage DS9
for images and TOPCAT for operations involving tables. But in astronomical
research, application software is usually not enough. For simple image
operations you might get by with firing up the DS9 software, for instance. But
at some point, sooner rather than later, you will want to do something more
specialised, and more automatised, than application software can provide.
Similarly, for data analysis. In some of the simplest cases, you might get
away with loading the catalogue in Microsoft Excel and start analyzing your
data in there. But in all other cases, including almost all of the interesting
ones, your analysis will need a little more flexibility. That is when, again,
you start writing a bit of code that helps you choose the right entries from
the catalogues, and to produce helpful diagrams – plots and histograms – that
allow you to make sense of your data. Then it becomes time to make use of a
different kind of tool:

  * •

A programming language is a tool for writing your own custom applications.

When you are using such a programming language for data analysis, you are in
effect writing yourself a custom application that can be used for the specific
analysis problem you need to solve. This approach has the advantage that you
have (nearly) full control over what you will be doing what your data.

It would not be an effective use of your time if you were to re-invent the
wheel by using a programming language for writing your own routines for
standard tasks such as opening different kinds of data files, or standard
analysis operations. Instead, you should make use of useful collections of
routine operations that have been written by helpful other people. In the
ecosystems of programming languages, these usually come in the form of
libraries or modules: chunks of codes that are easily included in your own
program, and give you pre-programmed functionality you can use for your own
specific purposes. That way, you need not write everything from scratch. But
you will still need to program in order to string these tools together to do
your bidding.

What qualifies as a routine operation will depend on context, of course.
Specialised astronomical modules provide you with tools for higher-level
operations that are routine in astronomy, but not elsewhere. An ephemeris
module will help you find the position of Solar System Bodies, for instance.
Some routines may be adapted to a specific telescope, allowing you to reduce
and analyze that telescope’s data. While you are still learning, you will want
to avoid some of those higer-level modules and re-invent at least some of the
wheels in question, since writing a routine for completing some specified
astronomical task is a good way of understanding what that particular task,
and the astronomy behind it. When have become more advanced, ready-made
modules represent a different problem: In research, it is important that you
understand the different steps in any analysis you are doing. Using a module
which is a “black box” for you represents a step in your analysis that you do
not fully understand. In those cases, it is even more important than usual for
your analysis to include suitable cross-checks and tests to ensure that it is
indeed doing what you intend it to do.

###  1.3 Concepts and operations

In general, when working with astronomical data analysis, you need not
understand all areas and aspects of a programming language. But there are
certain concepts, and certain types of operations/manipulations, which
constitute the basic working knowledge of virtually every astronomer working
with data.

This starts with basic mathematical operations. When you go from the magnitude
to the flux emitted by an astronomical object, you will need the “x to the
power of n” operation; on the way back, the logarithm. Whenever you perform
calculations with your data, you will need the appropriate operations.

Data points come in sets: the pixel data for an image corresponds to a two-
dimensional arrangement, while a list of properties for astronomical objects
will correspond to a one-dimensional chain of values. Programming languages
feature suitable data structures for this kind of connected data, such as
lists, arrays, tuples, or different kinds of table (the meaning of those words
can differ somewhat between one programming language and the next).

Knowledge of these data types and the various ways of manipulating them is a
must, along with knowledge of more basic types such as strings, integers or
floating point numbers — and of course the basic concept of storing values in
a variable in the first place!

For operations on our data, we need control structures. If we want to perform
a certain operation on every element of a list, for instance, we will need
something like a for loop. In order to distinguish between different cases — a
structure that allows us to apply a certain combination of operations to every
element of a list.

There will also be situations where we might want to perform a certain
operation on some elements of the list, but not on specific other elements —
to accomplish this, we need if clauses, that is, tools that tell our program
to apply certain operations only if specific conditions are met, but not
otherwise.

In addition to this kind of general knowledge, which is required when learning
pretty much any general programming language, astronomers should have at their
disposal a set of programming tools for more specific tasks — which often
equates with familiarity with particular libraries or modules. Often, we want
to visualize our data, so knowledge of how to create various kinds of plots,
diagrams or histograms (both one-dimensional histograms and two-dimensional
density plots) is essential.

Last but not least, how do we get data into our program, and our results out
again? If we have obtained the data by downloading a file, we will need to
know about proper input/output operations (in short, i/o). For certain data
formats, such as the ubiquitous FITS image files that are the usual format for
astronomical images, or for astronomical tables in FITS or VOTable format,
there are special functions that read the data in a way that makes it
particularly easy to start working with them.

When we do not download the data in the form of files, but instead access
astronomical data bases, there is an additional issue. We need to tell the
data base which specific set of data we would like to access. In order to do
so, we must submit a data base query, or query for short, to the data base: a
formalized request for data, written in a specific query language. A number of
astronomical data bases are organised in the shape of a Virtual Observatory
(VO) — data bases that conform to certain common standards to enable easy
access for all astronomers. The query language for the VO is the Astronomical
Data Query Language (ADQL), which is similar to the more general Structured
Query Language (SQL, pronounced either ”S–Q–L” or ”sequel”). Queries in this
language are useful both in the context of an application software like
TOPCAT, where they can used in the framework of the Table Access Protocol to
download a specific subset of data from an online data base via the Internet,
or as part of a Python program.

There is another aspect of all this, which would require a tutorial of its own
for proper treatment: data can be generated by software, too. Astronomy isn’t
only about observing. In the end, we want to understand the objects we
observe. That involves creating simplified models for these objects. If a star
is (put simply) a gigantic ball of plasma, held together by its own gravity
and heated up by nuclear fusion reactions in its core, then if we create a
simulation of such an object, using our knowledge of the laws of physics, the
resulting model should have similar properties to a real star (as we can check
using observations).

Simulations, too, require coding. In physics, only the simplest situations can
be described “analytically”, that is, writing down what happens in terms of
simple functions such as $\sin(x),\cos(x)$, polynomials and the like. For more
complicated situations, you will need to simulate what happens numerically:
starting with the initial situation, and then letting the computer
reconstruct, time step by time step, what happens. We will encounter a very
simple simulation in section 11.

###  1.4 Software/language choices

Every text on data processing has the same problem: For specific applications,
there is usually more than one application software, and of course there are
numerous programming languages. In teaching about data processing, one should
include specific examples, and students should work through such examples
themselves, gaining hands-on experience with all that data processing
involves. If the author chooses to present these examples in one specific
programming language, at least some students will later, when they are working
on a specific project, need to re-learn a different programming language.

This is not as bad as it sounds, though. Most programming languages, and most
applications, share similar concepts and functions. Once you have learned
about those in the framework of one specific programming language, or
application software, switching to another language or software will be much
easier than starting from scratch. Thus, learning what this tutorial has to
offer is definitely not a waste of time, even if it should turn out that later
on you will need to adapt to other software.

In this tutorial, I have chosen some common application software for simple
operations: SAOImage DS9 (DS9 for short) is a comparatively simple image
viewer that also allows some basic manipulation of astronomical images. TOPCAT
is a standard tool for accessing data from the Virtual Observatory. The
programming language used for more complex tasks is Python, which is widely
used in astronomy. This wide use has a great advantage: astronomers have been
writing helpful astronomy-specific libraries and modules for Python, and are
actively maintaining them. If you’re starting a career in astronomy, chances
are that you will do your basic programming in Python.

All that said, let’s get started. To get our bearings, we start with something
simple: before we delve into astronomical Python and start coding ourselves,
let us begin with two simple use cases for application software: In section 3,
we will look at astronomical images and combine red, green and blue filter
images into a color image. In section 4, we will look at some basic table
operations with TOPCAT.

##  2 Data basics: images, spectra, tables

In astronomy, just as in other sciences, we are not interested in data for
data’s sake. We want to do astrophysics: we want to use data to further our
understanding of the universe. Before we look at specific tools, and learn how
to use them, let us consider some of the properties of astronomical data, as
well as some of the specific ways of extracting information from them.

###  2.1 Images: Colour, brightness, pixels

Astronomers take images of astronomical objects, using telescopes and suitable
instruments attached to those telescopes. Public versions of such images can
be stunningly beautiful, and contribute significantly to the fascination of
the general public with astronomy. The underlying science images are commonly
stored in a format known as FITS, which stands for the “Flexible Image
Transport System” — a flexible file format that astronomers have been using
for images, spectra, data tables and more since the 1980s. When you encounter
a professional astronomical image, it is likely to be in that particular
format, with file extensions “.fits” or “.fit” on an older Windows machine.
Figure 1 shows one of the iconic images from the Hubble Space Telescope,
namely the open cluster Westerlund 2.

Figure 1: Image of the open star cluster Westerlund 2, taken by the Hubble
Space Telescope. The image data was downloaded from spacetelescope.org. Image
credit: NASA, ESA, the Hubble Heritage Team (STScI/AURA), A. Nota (ESA/STScI),
and the Westerlund 2 Science Team

Using this image as an example, we can demonstrate a number of properties of
astronomical image data. Phenomenologically, the image contains two types of
information: the stars we see in the image are much too small for even Hubble
to see any of their structure. They appear as point sources. In addition, we
have extended sources, in this case a region of ionized hydrogen (HII, in
astronomical parlance), which are, as the name says, extended areas with
varying brightness and colour.

While we tend to think of astronomical images as a rendition of “what’s up
there in the sky,” there are several aspects in which such images are not
faithful renditions — and those aspects are crucial for understanding
astronomical image data. Let us start with the colours. Professional
astronomical images are black-and-white images. One reason for this is that
digital cameras are, at their most basic, black-and-white. For each pixel,
they can only record how much light has fallen onto the collecting area for
that pixel (more specifically: how many photons have fallen). Consumer cameras
as in your smartphone or your digital camera are only able to produce colour
images because they have an array of filters installed in front of the array
of light-detecting sensor pixels. A common pattern is the Bayer mask, part of
which is shown in figure 2.

Figure 2: Part of a Bayer mask pattern: an array of filters installed in front
of detector pixels to enable the quick creation of a colour image

When such a consumer camera has taken an image, the colour information is
interpolated, and a colour image is displayed and saved. (If you have a camera
that can save images in some kind of “raw” mode, you can see the not-yet-
interpolated pattern.)

For astronomical images, a fixed filter mask is impractical for several
reasons. Astronomers would like to get the full resolution for their images,
so reducing resolution in each color band by having information about the
green brightness in every second pixel only, and about red and blue in every
fourth pixel, is a drawback. Colour interpolation means that some of the
colour information gets lost. Also, astronomers use a bewildering array of
possible filters, not just those corresponding to red, green, and blue (R, G,
B) — some of those filters capture a wider wavelength range, while narrow-band
filters might capture just a particular spectral line. Astronomers need the
flexibility of putting these different filters in front of their camera. On
the plus side, most astronomical objects change only very slowly. Using
different filters in succession, taking an image first through one filter,
then through the next, is perfectly feasible; those images will all show the
target object in effectively the same state. (On the rare occasions where
speed is of the essence, astronomers will use something like parallel cameras
observing through different filters. A case in point are gamma ray burst
afterglows, which fade on a time scale of several minutes. The GROND
instrument at the 2.2m ESO/MPG telescope at La Silla observatory was optimised
for observing such afterglows through seven different filters,
simultaneously.)

From these separate images, we can reconstruct colour images. (One simple way
of doing this will be shown in section 3.5.) If your aim is to produce a
“pretty picture,” there are several filter combinations that will serve.
Perhaps the best-known is the combination of images taken with the filters B,
V, and R respectively from the Johnson-Morgan photometric system to represent
the colours blue, green, and red.

But most images taken with the Hubble Space Telescope do not include B, V and
R versions. Colour images composed from them will be false-colour, and look
markedly different from the colours we see around us. That is OK as long as
you use the colours only to help you discern structural details, but you
should be aware that the result is not a faithful rendition of astronomical
colours. (There are heated discussions among amateur astronomers222These days,
with sophisticated technology available, the so-called amateur astronomers are
working very, very professionally indeed, and produce spectacular images.
about what constitutes faithful colours; some amateurs strive for the optimum
of faithful colour rendition. Their images are probably the closest you can
get these days to faithful colouring.) Figure 3 shows two images of the same
region side by side: the Westerlund 2 image already shown in figure 1, and an
image of the same region, taken in April 2017 with the 2 m Faulkes Telescope
operated by Las Cumbres Observatory at Siding Spring in Australia.

Figure 3: Two images of Westerlund 2. Image credit top: 2 m Faulkes Telescope
operated by Las Cumbres Observatory at Siding Spring under license CC BY-2.0.
Bottom image: NASA, ESA, the Hubble Heritage Team (STScI/AURA), A. Nota
(ESA/STScI), and the Westerlund 2 Science Team

The left-hand image was reconstructed from three images taken through B, V,
and R filters. The red colour of the gas cloud is probably fairly realistic.
The Hubble image no the right, while much higher resolved and full of
structural details, doesn’t show faithful colours. Few Hubble images do (which
does not take away from, but possibly adds to their power to fascinate
viewers).

Another aspect to keep in mind is the brightness scale of the image. FITS
images typically have 16-bit brightness values, that is, each pixel can take
on values from 0 to $2^{16}=65~{}536$ (apart from a possible offset value);
alternatively, brightness values within that range are stored as floating-
point numbers. Even in the integer case the best current computer screens
cannot faithfully display that dynamic range. And even if they could, the
result of simply displaying pixel brightness in a linear fashion would not let
you see the interesting details, and understand your image.

Published astronomical images will always333With possible exceptions when an
author is trying to make a point about scaling and astronomical images. have
some kind of (brightness) scaling applied to them. Figure 4 again shows the 2
m Faulkes Telescope image of Westerlund 2, this time only the one taken
through the red R filter. The leftmost version has linear scaling (pixel value
linearly related to brightness shown in the image), with the smallest pixel
value taken as black and the brightest taken as white. The only structure
visible corresponds to the locations of the brightest stars. No cloud
structure is visible in this rendition. The center version is still using a
linear map, but now every pixel value smaller than 4573 is displayed as black,
and every pixel value larger than 6001 is displayed as white, and brightness
values between 4572 and 6002 are displayed as the various shades of grey in
between.

Figure 4: Part of the Westlund 2 image taken through an R (red) filter in
April 2017 with the 2 m Faulkes Telescope operated by Las Cumbres Observatory
at Siding Spring in Australia with different scaling. Left: Linear scaling
from 0 to 65536. Center: Linear scaling from 4572 to 6002. Right: Square
scaling from 4572 to 6002

By concentrating on this narrow range of brightness values, we not only see
more stars, but also some of the structure of the cluster’s hydrogen clouds.
But at this scaling, the background sky looks rather bright. If we insert a
square function — where the displayed brightness is proportional to the square
of the pixel brightness value — we obtain a clearer differentiation between
the dark background and the brighter areas corresponding to the cloud. This is
shown in the version on the right.

Choosing a good scaling is not an exact science, but a matter of artisanship:
a good scaling will serve to illustrate the structures that are scientifically
interesting. But you should always keep in mind that there were choices
involved in creating the image.

A number of images are really mosaics, where several images have been stitched
together to form a larger picture. The beautiful Hubble version of Westerlund
2 in figure 1 is a case in point, as it is a composite image using
observations with Hubble’s Advanced Camera for Surveys (ACS) and its Wide
Field and Planetary Camera 3 (WFPC3). Figure 5 shows a sample WFPC3 image
(although probably not one used in the final composite444 I didn’t find any of
the original WFPC3 images in the MAST archive; all I could find were already
(smaller) composites.)

Figure 5: Image credit: NASA, ESA, the Hubble Heritage Team (STScI/AURA), A.
Nota (ESA/STScI), and the Westerlund 2 Science Team

pasted into the final colour image to give you an idea of the footprint of the
WFPC3. In fact, the WFPC3 inset is already a blend of four images, created
from the three image chips (CCDs) of the Wide Field Camera (the three larger
squares) and the image chip of the Planetary Camera (smaller square nestled
into the corner formed by the other three).

As a next step, let’s zoom in into the WFPC3 images, more concretely: into one
of the Wide Field Camera squares. The result can be seen in figure 6.

Figure 6: Zoom in on the WFPC2 image of Westerlund 2. Image credit: NASA, ESA

There are several points of note. The first is that the image is made up of
discrete square fields: pixels. That is no surprise if you have ever looked
very, very closely at digital photographs. It also means that, at the lowest
level, working with digital images means working with pixels, and with the
brightness value associated with each pixel.

Pixel positions are described by a pair of (integer) coordinates for each
pixel. A schematic example is shown in figure 7.

Figure 7: Pixels and pixel coordinates

The brightest pixel in the $5\times 5$ array would have the coordinates
$(4,3)$, since it is in the fourth column from the left, and in the third row
from the bottom. (Beware, other conventions exist! Some count rows from the
top. Some start the count at number 0, not number 1.)

###  2.2 Images: PSF and noise

Back to figure 6. The disk-shaped bright objects in the image are stars.
Here’s the thing: Westerlund 2 is at a distance of about 20 000 light-years
from us. At that distance, even an especially large super giant with 1500
solar radii should subtend an angle of a mere $0.002^{\prime\prime}$. Each
pixel in the WFPC2 image has a side length of $0.1^{\prime\prime}$. If our
image were a faithful map showing the exact direction whence light reaches us
from the sky, even these largest known stars would fall within a single pixel.
Instead, they and the much more common markedly smaller stars are smeared out
and appear in the image as disks. Figure 8 shows a brightness profile of the
star at the bottom center of figure 6. The disk is a few pixel wide.

Figure 8: Brightness profile of the bright star near the bottom center of
figure 6

Why the disk? Why not a single pixel? The answer, as you probably know, is
that light is a wave phenomenon, and that a wave passing through an opening —
in this case, the aperture of the telescope — is diffracted. The result is a
diffraction pattern that makes the image of a point source a disk (if you were
to look very closely, a disk with concentric rings around it). The larger the
telescope aperture, the smaller the disk — which is one key reason to build
ever larger telescopes: ever better resolution for the resulting images. The
function that defines the brightness distribution that results when a
telescope instrument produces an image of a point source is called the point-
spread function, abbreviated PSF.

Last but not least, look at the part of figure 6 that is not stars, but
background. The background is not uniformly black, but mottled grey — a
section, shown at even larger magnification, can be seen in figure 9.

Figure 9: Zooming in on part of the background of figure 6

There are several reasons that some of the pixels are brighter, others less
bright. One is the presence of distant, unresolved astronomical objects. But
those cannot explain the small-scale variation from pixel to pixel — remember
that even a point source would appear as a smeared-out disk! Instead, the
variability is noise — a spurious addition that tells us nothing about the
astronomical light sources out there.

The most fundamental effect is one of statistics: light reaches our detectors
in the form of photons, of light particles. The intensity of light reaching us
from a specific source determines the probability of a photon arriving within
a certain time interval. But the arrival itself is a random event. (You
probably know a similar situation: radioactive decay, where the decay
probability per unit time is constant, but each decay will still occur at an
unpredictable random time.) This randomness translates into pixel brightness
fluctuations. The relative size of these fluctuations shrinks as the total
number of photons collected grows. This is the other key reason why
astronomers want large telescopes (and in most cases still need long exposure
times, in addition): the more light they can collect from a distant object,
the smaller the relative fluctuations, the higher the signal-to-noise ratio,
and thus the clearer the image of those distant structures.

There are other kinds of noise. If you inspect raw, unprocessed images taken
with the Hubble Space Telescopes, you will find one kind that is typical for
space telescopes: traces left by cosmic particles depositing their energy in
the detector, leading to either longish streaks or more sharply defined dots,
depending on the direction the particle was travelling.

Figure 10: Traces of cosmic ray particles in an image taken with the Hubble
Space telescope. Image credit: NASA and ESA

Figure 10 shows an example (albeit from Hubble targeting a different object,
the Eagle nebula). Also, there is noise from the electronic devices involved
(although cooling key electronic elements down can reduce that kind of noise
considerably).

###  2.3 Images: Noise and flatfielding

When astronomers prepare data for the extraction of astronomical information,
a process commonly called data reduction, there are several typical steps they
take in order to reduce both the noise produced in their instrument and the
instrument’s idiosyncrasies when it comes to recording brightness.555In
preparing this section, I have profited from two online sources: The lecture
notes by S. Littlefair for the course PHY217 – Observational Techniques for
Astronomers he taught in 2014,
http://slittlefair.staff.shef.ac.uk/teaching/phy217/, last access 2019-11-01,
and parts of the e-book by R. A. Jansen, Astronomy with Charged [sic] Coupled
Devices (2006),
http://www.public.asu.edu/r̃jansen/ast598/2006ACCD.ebook…1J.pdf, last access
2019-11-01.

The image chips consist of little pixel elements; light falling onto a pixel
sets free some electrons. In a CMOS, each pixel also contains the electronics,
including a little amplifier, to read out a signal that indicates the number
of electrons, and thus the amount of light. In a CCD camera, the read-out
process is more involved, and involves herding electrons to the end of each
pixel row, then moving them to an amplifier. In both cases, ideally, the
number of electrons will be in direct proportion to the number of photons that
have fallen onto that pixel during the exposure time. And in the end, those
electrons are dumped onto a capacitor, whose voltage is measured. Since the
voltage across a capacitor is proportional to the accumulated charge, the
result gives us a measure of the number of electrons, and thus of the amount
of light we have captured. The analog voltage value is fed into an analog-to-
digital converter (ADC) which converts the voltage value into an integer
digital number, corresponding to analog-to-digital-units (ADUs) or, for short,
to counts associated with that pixel.

The conversion factor from the number of electrons to ADUs is called the gain,
typically given in units of $e^{-}$/ADU. Since the ADU values are digital, a
gain higher than 1 will introduce quantization noise through rounding errors —
for a gain of 6 $e^{-}$/ADU, the value 1 ADU could stand for anything between
1 and 5 electrons. A common range is for the digital numbers used to store the
count value to have 16 bits, and for each pixel to have count numbers between
0 and 65 535.

There is also the question of choosing the zero level for our voltage — that
is, the level that will correspond to the number 0 in our digitalized
documentation for each pixel. When no light falls onto the pixel, there is
still bound to be some level of noise, corresponding to a voltage value which
is called the bias. But mapping the average bias, call it $V_{0,avg}$, to the
zero value of the ADU is usually not a good idea. Our measurements of that
zero-light voltage are going to fluctuate slightly, its values sometimes a
little higher, sometimes a little lower. Frequently, we will want to average
images, either pixel by pixel or after aligning the astronomical objects in
those images. When taking such averages, above-average and below-average bias
levels should cancel out. But if we set zero ADU to correspond to the average
bias, then all fluctuations below that level will be mapped to zero as well,
and cannot cancel the fluctuations that are higher than average. You can take
care of the bias by taking a bias frame — an image taken with your camera
shutter closed (so no external light falls onto your chip) and an exposure
time as close to zero as your camera can manage. Subtract that bias frame —
or, again, the average of several such frames — from your science image, and
you’ve taken care of the bias.

Depending on the type of camera you are using, and on the type of
observations, you might need to perform different kinds of subtraction. For
older CCD cameras in particular, you might need to take a dark frame,
recording the contribution of electric currents due to thermal fluctuations,
and making a suitable correction. Under excellent conditions, and in
particular for infrared observations, you might need to take a sky frame to
subtract the sky brightness, including effects like airglow, faint radiation
emitted by the Earth’s atmosphere itself.

The corrections mentioned so far — bias frame, dark frame, sky frame — all
involved erroneous additions to pixel values, and the correction required the
subtraction of suitable correction frames. Another, different type of
correction is necessary because the sky is bound to get mapped to your digital
image unevenly — some parts of it too bright, others not quite bright enough.
I will illustrate this using a non-astronomical image, since the effect is
more obvious in an everyday setting. Fig. 11 shows a holiday snapshot, a lake
in Northern Germany.

Figure 11: Holiday snap with cxslearly visible vignetting

But you can see that the brightness is not as you would expect: near the
edges, the image becomes darker. This is not a property of Northern-German
lake landscape; it is a property of the camera I have used to record the
image, an effect known as vignetting.666Actually, I have added that effect
afterwards, artificially, but let’s pretend it is real. Real vignettes look
very similar. In a telescope image, you might also see dark, ring-shaped
smudges reminiscent of coffee stains; those are caused by dust flecks on the
telescope mirror.

Systematic brightness distortions like this can be removed as follows. In our
example in Fig. 11, the brightness of each pixel is a combination of the
brightness of the part of the object (in this case, landscape) we have
recorded on the one hand, and the sensitivity (or lack thereof) of our
telescope-instrument combination in that specific image region on the other.

We can reconstruct at least the relative sensitivity of our telescope-
instrument combination by taking an image of a scenery with completely uniform
brightness. In that case, all brightness variations are, by definition, not
due to our target scenery (which is completely uniform), but due to
sensitivity variations. Such an image can be seen in Fig. 12, and is called a
flat field (image). In astronomy, such images are produced either by pointing
the telescope at a uniformly lit canvas within the dome, or else by waiting
for dusk (or dawn) and pointing the telescope at the sky, whose scattered
light will be sufficiently uniform across the typical narrow fields of view of
astronomical telescopes (“sky flats,” specifically ”dusk flats” or ”dawn
flats”).

Figure 12: Flatfield image taken with the same telescope-instrument
combination as the holiday snap, with clearly visible vignetting

For the flatfield image, too, one can produce dark frames (taken at the same
exposure as the flatfield image) and subtract them, resulting in a suitably
corrected master flat. That master flat encodes the sensitivity for each image
pixel. With this information, we can correct for the sensitivity variations as
follows. Assume that a pixel in the master flat is twice as bright as a second
pixel. That would mean our setup is only half as sensitive for the second
pixel than for the first. But if we were to take an image, and then to
multiply the brightness value for the second pixel with the factor two, we
would have compensated for the different sensitivity levels.

More generally, we can restore proper relative brightnesses of all our pixels
by dividing our science frame by our master flat, pixel by pixel. The result
for our holiday snap can be seen in Fig. 13.

Figure 13: Flatfield-corrected holiday snap

The compensation is not perfect. For instance, if fewer photons have reached a
certain detector pixel, then the statistical noise will be somewhat larger for
that pixel. Dividing the pixel value by a factor, as one does in flatfielding,
will not get rid of that additional noise.

All in all, we have learned the basics of how astronomers are reducing their
image data — compensating for noise that is added to each pixel by subtracting
suitable compensation terms (bias frame, dark frame, sky frame), and
afterwards compensating for sensitivity variations by dividing by a suitable
compensation term (flatfield image).

To sum up the last few sections: The digital astronomical images used by
astronomers are made of pixels; what we see is in part determined by the
properties of our target object, but in part by the properties of the optical
system used (telescope plus instrument and their PSF), and in part by noise.
The “elementary images” are black-and-white, and usually taken through a
specific filter. When such images are displayed, additional decisions were
involved about how to represent brightness. Published images frequently use
colour to convey additional information — although in most cases, these are
false colour images, which do not reproduce the colour of the object we would
perceive could we view it directly. Astronomers employ dark frames and
flatfielding to reduce certain types of noise, and of sensitivity variations.
Sometimes, images are fit together to form a larger mosaic.

When working with images, we need to keep all this in mind — after all, we
want to use the information contained in the image to make deductions about
the astronomical objects observed. To do that, we need to know which aspects
of the image really do contain information about the object — and not
information about the telescope-instrument combination, or photon statistics.

###  2.4 Images: astronomical information

So far, we have talked mostly about image artefacts — what makes an image
different from the real thing. Time to talk about the physics behind it all:
What information is contained in astronomical images? The information
important for classical astronomy, for a start: Images contain position
information about astronomical objects, information about where exactly an
object is located in the sky (for object whose position does not change in the
usual coordinate system), or about how its position changes over time.

In the era of classical astronomy, this was the main purpose of observatories:
determining the positions of stars in the sky, as an aid to celestial
navigation. This also included measurements that allowed for precise time-
keeping: until the advent of stable quartz clocks in the mid-20th century,
documenting the periodic changes in the night sky, in particular the diurnal
motion during one (sidereal) day, was the most accurate time-keeping method.

In modern astronomy, determining stellar positions remains an important sub-
field, which for the last few years has been dominated by ESA’s astrometry
satellite Gaia. Accurate catalogues of stellar positions not only provide a
framework for localising astronomical objects in general. Via the parallax
effect, they also provide information about the distances of stars in our
cosmic neighbourhood, which in turn is a prerequisite for farther-reaching
methods of astronomical distance determination. Knowledge of astronomical
distances is crucial for making deductions about object’s luminosity. (In
principle, an object that appears to be bright in the night sky could have a
rather faint luminosity, but appear bright since it is very close to us, or
else have a really high luminosity while being rather more distant.)

Then, there is photometry, that is, determining the (apparent) brightness of
astronomical objects. If we have chosen proper exposure times, and made all
necessary corrections, we could deduce the brightness of a star from the sum
of the values of the pixels associated with that star. In practice, it’s
difficult to separate star pixels from non-star pixels, but there is a simpler
way known as aperture photometry:

Figure 14: Aperture photometry in a part of the Hubble Space Telescope image
of Westerlund 2. Image credit: NASA, ESA, the Hubble Heritage Team
(STScI/AURA), A. Nota (ESA/STScI), and the Westerlund 2 Science Team

Define a circular region that contains the PSF of the star completely (the
yellow circle in Fig. 14). At some small distance outside, define an annular
region (bounded by the two blue circles in Fig. 14). Assuming that the central
circle contains only the star we are interested in, and the annular region no
discernible star at all, we can argue as follows: If we sum up the pixel
values within our central circle, we get the light from the star plus
background light. We can estimate the background light as follows: On average,
the background brightness should be the same in the central circle and in the
surrounding annulus. In the annulus, we can determine the average brightness
by summing up the pixel values in the annulus and then dividing by the number
of pixels (equivalently, by the area of the annulus). Multiply this average by
the number of pixels in the central circle (equivalently, by the area of the
circle), and the result will be the brightness contribution from the
background within the central circle. Subtract this contribution from the
total sum of the pixel values within the central circle, and what is left is a
measure of the brightness of the star. Note that photometric measurements are
sometimes made with the telescope slightly out of focus, distributing the
object’s light over a greater number of detector pixels for greater accuracy.

In astronomical practice, stars are point-like objects. For extended objects,
we can measure a surface brightness, given in brightness per angular area in
the sky, and we can measure how that brightness varies from location to
location. Such brightness maps contain information about the amount of
material we are seeing. The situation is more complicated when densities are
so high that some of the matter obscures our view of whatever matter lies
behind (that is, if the matter in question is “optically thick”). In the
simplest case, we can see all of the light from the matter of, say, a nebula
(the nebula is “optically thin”), and the brightness in a certain area of the
sky allows us to estimate how many atoms we are seeing in that area — a column
density since we cannot deduce the three-dimensional structure, only the
number of atoms within that column of the three-dimensional object which gets
projected to the sky-region in question.

Brightness measurements will only ever cover some limited region of the
electromagnetic spectrum. Some of the limitation comes about by the kind of
telescope we use. An ordinary optical telescope will be able to receive
visible light, near-infrared light, and ultraviolet light (which, for ground-
based telescopes, is somewhat pointless since almost all UV light is filtered
out in the Earth’s atmosphere). But its camera would not be able to detect,
say, mid-infrared light, let alone X-rays. In practice, as we have already
seen, astronomers voluntarily restrict themselves to even narrower portions of
the spectrum, by using suitable filters. This allows for quantitative
description of the colors of astronomical objects. An object that is bright
when viewed through a blue filter, but dim when viewed through a red filter,
will be blueish in color.

What we have called the brightness so far, summing up pixel values in our
image, is proportional to the number of photons from a certain source (or a
certain area of the sky) entering our telescope during the exposure time.
Since the exposure time is the same for all the objects in our image, the
ratio of brightness values for two such objects is equal to the ratio of the
energy per unit time (in the given filter band) we receive from those two
objects.

Also, since in both case we are using the same telescope, and hence the same
collecting area, the ratio is equal to the ratio of the energy per unit time
per unit area (again in the given filter band) for those objects, or using the
appropriate technical term: the ratio of their fluxes.

Flux ratios are how astronomers traditionally compare the apparent brightness
of celestial objects — except that, to ensure some degree of backwards-
compatibility with the naked-eye-based, 2000-year-old Ancient Greek magnitude
system (as one does), those ratios are measured on a logarithmic scale.
Specifically, if $F_{1}$ and $F_{2}$ are the flux values for light we receive
from two objects 1 and 2 in a specific filter band, then their apparent
magnitudes in that band are defined as

$m_{1}-m_{2}=-2.5\cdot\log\left(\frac{F_{1}}{F_{2}}\right).$ (1)

A reference point for the magnitude system is chosen by setting a value for
the magnitude of a specific star in a specific filter band; for instance, in
the V filter band that roughly corresponds to a green filter, the star Vega
was originally chosen as a zero point, although his modern visual magnitude is
$m_{V}=+0.03$.

Note the minus sign — magnitude values are larger for fainter stars. With the
naked eye, under good conditions, you can observe stars with magnitudes up to
about $m_{V}=6.5.$ For extended objects, we would need to document the energy
per unit time and unit collecting area that reaches us from a given solid
angle in the sky. There, the intensity, as energy received per unit time per
unit area per unit solid angle is the appropriate descriptive quantity.

Brightness, of course, can change. Different types of variable stars, for
instance, can be distinguished by the shapes of their light curves, which
document how their brightness changes over time. The transit method for
detecting exoplanets also relies on light-curve measurements.

Last but not least, images also contain spatial, 3D information about the
objects under study. Typically, that information is projected onto the sky —
we do not see the full three-dimensional structure of, say, a gas cloud;
instead, we have one particular fixed perspective on that cloud. Interpreting
what we see typically involves models for the physical, three-dimensional
structure, whose predictions can then be compared to what we actually observe.

###  2.5 Spectra

On to a central kind of data set in astrophysics: spectra! A spectrograph
contains a dispersive element (or even more than one), which splits the
incoming light into its rainbow colours or, in physics terminology, into its
different wavelengths. Examples for the three basic types of dispersive
element an be seen in Fig. 15.

Figure 15: The three basic types of dispersive element, from left to right:
transmission grating, prism, and reflective grating

Figure 16 shows ceiling lamps, imaged through a dispersive grid, namely
through “spectral glasses” that can be used to demonstrate dispersion effects.

Figure 16: Light sources, imaged through a dispersive element (“spectral
glasses”)

As you can see in the figure, though, the spectral decomposition makes for a
hodgepodge of effects. The coloured images of the lamps, each corresponding to
a spectral line, overlap, creating a mix of spatial information and wavelength
information that is not easily disentangled.

Figure 17: One of the light sources masked by a slit, and imaged through the
same dispersive element

A common solution is shown in figure 17: introduce a slit mask, with the slit
oriented at right angles to the direction of spectral dispersion. If you keep
the slit small, the overlaps between the successive images of the slit will be
small as well. In the case of separate spectral emission lines, as in this
example, you can then see the lines clearly separated.

When you see the image of a spectrum that is a broad colourful band (possibly
with bright emission lines, or dark absorption lines), you see a succession of
such slit images, each copy indicating the brightness at that particular
wavelength. One such image can be seen in figure 18: a spectrum of sunlight
reflected by the Moon, produced with a Baader DADOS spectrograph.

Figure 18: Spectrum of sunlight, reflected by the Moon

The image was taken for demonstration purposes with a consumer digital camera,
and thus shows colours — any professional image of a spectrum would, of
course, be in black and white; the colour carries no additional information,
as the position of a spectral feature along the horizontal axis already
defines its wavelength, and hence its colour.

The information contained in a spectrum is one-dimensional — for each
wavelength, we have a quantity indicating how much light is emitted in that
particular wavelength region. Thus, it is natural to plot spectral data as a
curve. The top part of figure 19 shows an artificial color image of a Solar
spectrum, complete with dark Fraunhofer lines. The bottom plots the spectrum
as a curve, with wavelength plotted along the x axis and relative flux on the
y axis. The requisite data is taken from IAG solar flux atlas777 Reiners et
al. 2016,
[http://adsabs.harvard.edu/abs/2016A&A…587A..65R].

Figure 19: Solar spectrum as a curve (bottom) and the reconstructed version of
a slit spectrum image (top). Data from the IAG Solar Flux Atlas, Reiners et
al. 2016

The curve is quite complex, with a forest of absorption lines — narrowly
defined minima — one next to the other. If we zoom in by plotting a much
smaller wavelength interval, as in figure 20, you can see the the lines
themselves have characteristic shapes.

Figure 20: Narrow region from a solar spectrum. Data from the IAG Solar Flux
Atlas, Reiners et al. 2016

When working with raw data from spectra, certain reduction steps need to be
taken. Some of those are similar to the reduction of image data: Flatfielding
is again needed to compensate for differing sensitivity of the instrument in
different parts of the spectrum. This is more difficult for a spectrum than
for an image, since for a true flat field, you would need a perfectly flat
spectrum. Instead, any well-known, preferably smooth calibration spectrum can
be used to deduce the varying sensitivity. Dark frames again can be used to
take into account that the electronics of the detector will produce some
spurious brightness in the image, which needs to be subtracted.

Wavelength calibration is another necessity. After all, the spectral spread
has a specific meaning — light is separated according to wavelength (or
frequency). In order to map specific wavelengths to the direction along which
dispersion takes place, astronomers often employ specific calibration lamps,
which contain a gas or a mixture of gases that produce a hopefully dense array
of known emission lines. For a simple amateur spectrograph, you might use Neon
for the purpose; at the professional level, you might for instance find a
mixture of Thorium, Argon, and Neon. Sometimes, the calibration lines will be
recorded separately; in other cases, they are recorded concurrently with the
astronomical observation to allow for a direct comparison. A special case of
the latter, unavoidable for ground-based telescopes, are telluric lines —
absorption lines created not in outer space, but by light absorption in the
Earth’s atmosphere. Such telluric lines can be used for calibration, as well.

Image distortions can make spectral data reduction particularly challenging. A
particularly complex case are Echelle spectrographs, where two kinds of
spectral dispersion are combined: A grating will, in fact, produce several
different spectra, called “spectra of different order”. Higher orders tend to
overlap each other, but the different overlapping partial spectra can be
separated by a second stage of dispersion in the direction orthogonal to the
initial dispersion. The result are different rows of partial spectra, allowing
astronomers to capture an entire high-resolution spectrum on a standard, two-
dimensional camera chip. The raw image of one such Echelle spectrum, taken
with the FEROS spectrograph at the MPG/ESO 2.2-metre telescope at ESO’s La
Silla observatory, can be seen in Fig. 21.

Figure 21: Region within a raw image of a spectrum of the star HIP66974 and a
calibration lamp, taken with ESO’s FEROS spectrograph in June 2015 (Data set
SAF+FEROS.2015-06-13T23:16:46.772). Retrieved from ESO’s Science Archive on 24
April 2019

This is just a small region, about 10%, of a much larger image. As you can
see, the horizontal, curved stripes always come in pairs: the stripe on top is
mostly white with some dark absorption lines (which increase lower in the
image), while the lower stripe consists of a fairly dense forest of emission
lines. The upper stripe of each pair is the science image, in this case of the
star HIP66974, a star with the same spectral type (G2V) as the Sun and thus a
fairly similar spectrum.

The lower stripe in each instance is the calibration lamp — hence the many
emission lines, each marking a well-known reference wavelength. Reducing this
spectrum would mean to map the different stripes to their proper wavelength
regions (using the calibration lines), unbending the curved image, and
properly calibrating the brightness over the different parts of the image.

Once the spectrum is reduced, or if one is working with a reduced spectrum in
the first place, the spectrum as a whole and in particular the spectral lines
contain a wealth of information about the object in question. Systematic
Doppler shifts in the spectral lines indicate whether or not the light source
is moving towards us or away from us. Doppler shifts that change periodically
over time contain information about objects orbiting each other, from double
stars to exoplanets detected by the radial velocity method. Simple data
analysis in these cases proceeds by fitting the individual spectral lines,
finding their central wavelength, and tracing the changes of that wavelength
over time.

The shape and relative depth of spectral lines of a star contains information
about the star’s metallicity, that is, the fraction of elements heavier than
Helium contained in the star’s atmosphere, about the surface gravity and about
the effective temperature. Some specific spectral lines corresponding to
radioactive elements can be used to reconstruct the age of stars, and have
been used to find the oldest stars in existence.

The simplest part of such an analysis is about identifying the lines
corresponding to specific chemical elements; these lines show which elements
are present in the star’s atmosphere.

Some such lines are indicated in figure 19. In all of these cases, analysis
usually proceeds by creating spectra based on suitable models and comparing
those with the actual observations, finding the best fit.

###  2.6 Data cubes

So far, we have talked about two-dimensional images (where the two dimensions
correspond to an area on the night sky) and one-dimensional spectra (where the
one dimension corresponds to wavelengths). Data cubes are the combination of
this: We have a two-dimensional image of the night sky, but at each pixel
location, we have not only a single brightness value, but instead a whole
spectrum.

With two plus one dimensions, we are effectively looking not at a two-
dimensional rectangle, but a three-dimensional cube. A data cube does not
contain all the information reaching us from a certain region of the sky at a
certain time (polarisation information is missing), but it comes impressively
close.

One way of obtaining such a data cube is with Integral Field Spectroscopy —
for instance: splitting an image into comparatively large “pixels,” each of
which is channeled into a glas fibre which transmits its light to a
spectrograph, where the spectrum is then recorded.

Another natural way of recording such a data cube is in interferometry, where
the light from several telescopes is combined in a coherent way, making use of
the wave nature of light. In reconstructing images from interferometric
measurements, one can distinguish (to a certain degree) between contributions
with different wavelengths; in effect, this allows for the reconstruction of a
three-dimensional data cube.

Human beings are not equipped for really three-dimensional vision. What we
call three-dimensional vision is really just seeing surfaces within (sparsely
populated) three-dimensional space. We cannot see all the points within a
three-dimensional data cube at once.

Figure 22: 16 of the 72 channels recorded for the galaxy NGC 3198, from the
THINGS survey (Walter et al. 2008), [http://www.mpia.de/THINGS/]

Fig. 22 shows one solution: showing separate images for different regions
within the spectrum. In this particular case, the “channel map” shows 16 of
the 72 frequency bins around the 21 cm hydrogen line that is characteristic
for atomic hydrogen (that is, hydrogen atoms; not bound into hydrogen
molecules, not ionized to form a plasma, just simple atoms).

One of the most interesting applications of data cube data is to extract the
information they contain about the large-scale motion of matter. Fig. 22 shows
21 cm radiation emitted by hydrogen atoms, but some of that radiation is
shifted to lower and some to higher frequencies — why? Because some of the
atoms are moving towards us, others away from us, and their 21 cm radiation
undergoes a corresponding Doppler shift.

The data cube contains information on the radial velocity of the gas we see in
the different frequency channels. We can combine that information to make a
color picture whose color encodes the average radial velocity of gas in each
region of the image, giving what is called a first moment map. That picture is
shown in Fig. 23.

Figure 23: Velocity map for the galaxy NGC 3198, from the THINGS survey
(Walter et al. 2008), [http://www.mpia.de/THINGS/]

Reddish regions are moving away from us, blueish regions towards us. The
combined picture is that of a disk galaxy whose stars are rotating as a whole,
one side of the disk coming towards us and the other moving away from us in
that coordinated motion.

Alternatively, instead of taking the average, we could compute the standard
deviation of the velocity values associated with each pixel. That would give
us an estimate not of the bulge motion of gas in that region, but of the
diversity of motion, the spread of radial velocities.

In a similar manner, we can use data cube information to map all those
quantities that can be derived from spectra — the presence of specific
elements and thus the chemical composition, on larger scales the prevalence of
different kinds of stars, and more.

###  2.7 High-level data: catalogues and tables

Once astronomers have derived observational or physical quantities from their
observations and measurements — deduced the temperature of stars from their
spectra, or their luminosity from their apparent brightness and some measure
of their distance — they can compile catalogues containing such higher-level,
derived physical information. Analysing this kind of high-level data is
broadly similar to statistical analysis in other fields, and uses the same
general tools.

To begin with, a catalogue is no more than a list. Conventionally, each row in
that list represents a separate object, and each column represents a property.
If you just look at the numbers in a big list, you are sure to miss the forest
for the trees. In order to extract trends, distribution, systematic
correlations, functional relationships from the data, in astronomy: in order
to understand what physical laws and evolutionary pathways has made objects
the way they are, we need employ proper tools. Statistical analysis is a wide
field, and in this basic introduction, we will only look at the most basic of
descriptive tools.

A very basic tool is a histogram, which allows us to see the distribution of
values for a single quantity within our sample (e.g. among the objects of our
catalogue). The basic principle is simple: within the range $x_{min}\leq x\leq
x_{max}$ spanned by the values for the quantity $x$, we define $N$ bins of
equal size, each with a lower boundary $x_{i}$ and an upper boundary
$x_{i+1}$. A value $x$ falls into the bin with index $i$ if $x_{i}<x\leq
x_{i+1}$. We then draw the bins as rectangles, whose height is proportional to
the number of values which fell into the bin. This gives us a measure of how
prevalent (or not) specific values are.

Histograms can also be logarithmic. In that case, we divide the range of the
quantity we intend to map into bins of equal logarithmic size. For instance,
if one bin contains stars that are between $1=10^{0}$ times and $10=10^{1}$
times as luminous as the Sun, the next bin would be from $10^{1}$ to $10^{2}$
solar luminosities $L_{\odot}$, and the following one from $10^{2}$ to
$10^{3}\>L_{\odot}$. Those bins are of equal size when it comes to their
exponents. Such logarithmic binning is useful for physical properties that are
spread across a wide spectrum of scales.

Here is a simple example for a histogram with logarithmic bins. Fig. 24 shows
the basic preparations for manually drawing a histogram: a tally sheet for
putting the stars in a certain data set into their proper bins. The visual
appearance of the tally sheet already constitutes a simple histogram-like
representation, although in an unusual sideways orientation.

$10^{-3}$ to | $3\cdot 10^{-3}\;L_{\odot}$ |  | 0
---|---|---|---
$3\cdot 10^{-3}$ to | $10^{-2}\;L_{\odot}$ | 𝍷𝍷𝍷 | 3
$10^{-2}$ to | $3\cdot 10^{-2}\;L_{\odot}$ | 𝍸𝍸 | 10
$3\cdot 10^{-2}$ to | $10^{-1}\;L_{\odot}$ | 𝍸𝍸𝍷𝍷𝍷 | 13
$10^{-1}$ to | $3\cdot 10^{-1}\;L_{\odot}$ | 𝍸𝍸𝍷𝍷 | 12
$3\cdot 10^{-1}$ to | $10^{0\phantom{-}}\;L_{\odot}$ | 𝍸𝍸𝍸𝍸𝍸𝍷 | 26
$10^{0\phantom{-}}$ to | $3\cdot 10^{0\phantom{-}}\;L_{\odot}$ | 𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍸 | 60
$3\cdot 10^{0\phantom{-}}$ to | $10^{1\phantom{-}}\;L_{\odot}$ | 𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍸 | 55
$10^{1\phantom{-}}$ to | $3\cdot 10^{1\phantom{-}}\;L_{\odot}$ | 𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍸𝍷𝍷𝍷𝍷 | 44
$3\cdot 10^{1\phantom{-}}$ to | $10^{2\phantom{-}}\;L_{\odot}$ | 𝍸𝍸𝍸𝍸𝍸𝍸𝍷𝍷𝍷𝍷 | 34
$10^{2\phantom{-}}$ to | $3\cdot 10^{2\phantom{-}}\;L_{\odot}$ | 𝍸𝍸 | 10
$3\cdot 10^{2\phantom{-}}$ to | $10^{3\phantom{-}}\;L_{\odot}$ | 𝍸𝍸𝍸𝍷𝍷𝍷 | 18
$10^{3\phantom{-}}$ to | $3\cdot 10^{3\phantom{-}}\;L_{\odot}$ | 𝍸𝍸𝍷 | 11
$3\cdot 10^{3\phantom{-}}$ to | $10^{4\phantom{-}}\;L_{\odot}$ | 𝍸𝍷 | 6
$10^{4\phantom{-}}$ to | $3\cdot 10^{4\phantom{-}}\;L_{\odot}$ | 𝍸 | 5
$3\cdot 10^{4\phantom{-}}$ to | $10^{5\phantom{-}}\;L_{\odot}$ | 𝍸𝍷 | 6
$10^{5\phantom{-}}$ to | $3\cdot 10^{5\phantom{-}}\;L_{\odot}$ | 𝍷𝍷𝍷 | 3
$3\cdot 10^{5\phantom{-}}$ to | $10^{6\phantom{-}}\;L_{\odot}$ |  | 0
Figure 24: Tally sheet in preparation for a hand-crafted histogram

Hand-crafted histograms have gone the way of so many other hand-crafted
things. The automatically-generated version for the same histogram can be seen
in Fig. 25.

Figure 25: Automatically-generated histogram

We will look at ways of generating such histograms in sections 4.7 (TOPCAT)
and 8.7 (Python).

The data set is DEBCat,888Southworth 2014, [https://arxiv.org/abs/1411.1219].
The data and supplemental information are available online at
[http://www.astro.keele.ac.uk/jkt/debcat/]. a collection of more than a
hundred well-studied transiting double stars. In that particular set-up,
astronomers have sufficient information to be able to reconstruct the stars
key physical properties like mass, radius, luminosity and temperature, and
derive properties such as a star’s average density.

Let us use some of those properties to look at another very common tools are
diagrams populated with data points in which we plot one quantity against
another. Each axis stands for a specific physical quantity (commonly scaled
either linearly or logarithmically), and each data point corresponds to a pair
of values, one for each of the axis quantities.

A famous example, with two logarithmically scaled axes, is the Hertzsprung-
Russell diagram, a version of which is shown in Fig. 26.

Figure 26: Hertzsprung-Russell diagram using physical quantities (temperature
and luminosity) instead of spectral classes and luminosity classes. Data from
DEBCat

The values corresponding to each data point are read off in the usual way, as
shown by the two grey auxiliary lines: from the dot, go horizontally to the
vertical axis; read off the value indicated at the intersection point, in this
case a luminosity of about $5\cdot 10^{4}\>L_{\odot}$ and a temperature of
somewhat more than $20~{}000\>\mbox{K}$.

In the diagram, you can clearly discern a linear structure going from bottom
left to top right, and a cloud of dots hovering to the upper left of that
linear structure. Astronomers call the linear structure the main sequence, and
the stars in it the main sequence stars.

There is a way of include additional information in histograms and
2-dimensional diagrams: use colour! In the simplest case, we can use colour to
distinguish between different populations of data points. For instance, let us
colour the points in the upper-left cloud on the Hertzsprung-Russell diagram
red, as Fig. 27.

Figure 27: Hertzsprung-Russell diagram using physical quantities (temperature
and luminosity) instead of spectral classes and luminosity classes. Data from
DEBCat

So far, the colouring hasn’t brought us any great advantage. The cloud was
apart from the rest before, and it is apart from the rest now. But let us
carry this color scheme over to a histogram, for instance, plotting histograms
for the red and the blue dots side by side, using the same bins. The result
for a histogram of radii is shown in Fig. 28.

Figure 28: Separate histograms for the blue and red data points from Fig.
HRDiagramRedBlue

Now we see that the stars corresponding to those red dots have considerably
larger radii than their main sequence counter parts. The size distributions
are clearly separate. Those red-dot stars are veritable giants! We know from
the Hertzsprung-Russell diagram that their temperatures are somewhere between
4000 and 6000 K, going from reddish to yellowish. So these red giant stars
were named with excellent reason.

Let’s look at a diagram plotting, say, radius against density, as in Fig. 29.

Figure 29: Plotting radius against density, both in solar units. Data from
DEBCat

In that diagram, it is not clear which are the main sequence stars and which
are the red giants. With the red-blue distinction, the situation becomes
clear, as shown in Fig. 30.

Figure 30: Plotting radius against density, both in solar units. Color marks
red giants vs. main sequence stars. Data from DEBCat

Red giants are not only generally larger in radius than main sequence stars,
they are also considerably less dense. Our data points us in the right
direction: in the modern view, main sequence stars go through a red giant
phase after they exhausted the hydrogen fusion fuel in their cores, their
atmospheres swelling up and cooling down in the process, leading to a large,
reddish star with drastically reduced mean density.

We can use colour more quantitatively than just to express class membership in
a two-colour scheme. Colour can add an (imperfect) third dimension to our
diagrams. In the version of a mass-luminosity diagram shown in Fig. 31, each
data point has the proper star color corresponding to it’s temperature (as
determined from the star’s spectral properties999Information about this kind
of color mapping can be found on
http://www.vendian.org/mncharity/dir3/starcolor/). This color-coding
immediately allows you to identify the red giants, see that their mass range
is a subset of the mass range of the main sequence stars, but that the red
giants are larger and more reddish.

Figure 31: Mass-luminosity diagram with data points plotted in the color
corresponding to a star’s temperature. Data from DEBCat

Color scales can also be artificial, and different color maps are available
for the purpose. In Fig. 32, the color now indicates the radius of each star,
with the scale shown by the colorbar on the right.

Figure 32: Mass-luminosity diagram with data points plotted using an
artificial color map that indicates the stars’ radii. Data from DEBCat

Clearly, stellar radii grow along the main sequence, but the red giants, in
their little cloud of data points above the main sequence, are larger still.

A good color map can make your diagram much easier to understand; a bad one
can be confusing. Also, you should take into account accessibility issues.
Your color maps should be accessible even to people with certain forms of
colour-blindness.101010Some information on this kind of accessibility can be
found on [https://betterfigures.org/2015/06/23/picking-a-colour-scale-for-
scientific-graphics/] Unless you are plotting a spectrum, avoid the rainbow
colour map. Instead, consider color maps that have been designed to be
accessible for those with colour blindness, as well as to print well in black
and white — for instance the Viridis family of colour maps.

Patterns in such diagrams indicate interesting relationships. Is there a
linear relation — do data points for certain physical quantities fall on a
straight line in a linear diagram? Or is there a power law at work, $y\sim
x^{a},$ in which case the data points would fall on a straight line in a log-
log diagram? In this way, diagrams can help us find systematic relations
between our data.

This is not as straightforward as it sounds, of course. In a two-dimensional
diagram, we can plot at most three different quantities (if we make clever use
of a color map). We could try all different pairs of physical quantities
relevant for the situation we are looking at, and might get lucky in finding
interesting relationships in that way. A three-dimensional diagram is
possible, but would need to be interactive so we can view it from all
different sides to get a feeling for the three-dimensional structures. Of
course, the basic physical quantities can be combined to yield compound
quantities. Complex relationships between quantities, longer polynomials
involving several quantities for instance, or differential/integral relations,
are much less straightforward to read off such diagrams.

The typical way of extracting information about systematic relationships from
data is to fit a function to the data. Assume that we have data points
$(x_{i},y_{i})$ for $i=1,\ldots,N$, each representing a pair of quantities. A
common measure for how well those data points satisfy a general relationship
$y=f(x)$ is as follows.

If the relationship were to hold perfectly, then we would have
$y_{i}=f(x_{i})$. In real life, functional relationships are not that perfect.
Even in cases where the relationship $y=f(x)$ is the basis for our set of data
points, measurement errors will lead to deviations. Moreover, exact
relationships are rare; the much more common case is that the relationship is
approximate, and that data points scatter around the curve $y=f(x)$.

For a single data point, the quantity

$\Delta y_{i}=y_{i}-f(x_{i})$ (2)

is a measure of the deviation of the data point from the relationship. What is
the best way of summarising the deviations for our data set as a whole? We can
say what is definitely not a good measure: taking the sum of all the $\Delta
y_{i}$, since deviations may be positive and negative, and the sizeable
deviations associated with different data points could cancel each other out,
skewing the result — we could even get an overall measure of zero, indicating
no deviation, in a situation where the $\Delta y_{i}$ are huge, but cancel
pair-wise!

To avoid this, we could take the sum of the absolute values $|\Delta y_{i}|$
but as we shall see later on, it is useful for the measure we choose to be
differentiable. That why a better choice is the sum over the squares of the
deviations,

$S\equiv\sum_{i=1}^{N}[y_{i}-f(x_{i})]^{2}.$ (3)

Commonly, we have an idea for the basic properties of the function $f(x)$, but
not about the explicit form of the function. For instance, we might have
reason to believe the function to be linear (since that is what an x-y diagram
suggests), $f(x)=ax+b,$ but do not know the values for $a$ and $b$.

In such a case, we can use the quantity (3) to find the best fit. Let us make
explicit that the function depends on the parameters $a,b$, namely as
$f(x,a,b)=ax+b.$ For our set of data points and for any given pair of values
$a,b$, we have

$S(a,b)\equiv\sum_{i=1}^{N}[y_{i}-f(x_{i},a,b)]^{2}.$ (4)

Fitting the function to the data involves minimising $S(a,b)$; since the
expression in question is the sum of the squares of the deviations, this is
known as the least-squares method. For linear functions, there is even an
analytic solution: At the minimum of $S(a,b)$ as a function of $a$, the
derivative of $S(a,b)$ with respect to $a$ must be zero (this is why it was
useful to choose $S(a,b)$ to be readily differentiable). The same goes for
$b$; impose both of those conditions, and you can find $a$ and $b$ directly.
In the general case, no such analytical solutions are possible, and the best
fit is found numerically.

There are a number of ways to go from here, some more advanced, some less. In
fitting a function to data points, you can give different weight to different
contributions to (4) with the measurement errors for each data point; in this
way, those data points that are less-well known will also contribute less to
the choice of parameters.

Then, there is the problem of outliers, that is, lonely data points that are
far from the rest, and probably not because of their physical properties but
because of measurement errors. Astronomers also make use of Bayesian
techniques in order to estimate the parameter values best fitted to their
model, or to compare different models with each other. All of these issues are
beyond the scope of at least this version of my basic introduction.111111Some
additional information can be found in the classic paper by Hogg, Bovy and
Lang 2010, [https://arxiv.org/abs/1008.4686].

The problem of making sense of multidimensional data is of interest far beyond
astronomy, and a main task of what has become known as data science.
Astronomers apply numerous tools that have much more general application in
order to solve their data science problems.

For astronomers, this is not a matter of merely applying well-established
methods and tools to new data sets. Instead, the way that astronomers look at
their data, and make their deductions, is continually evolving. It is
certainly more common for astronomers to derive new results by utilizing new
data, but it is equally possible to make new deductions from an existing data
set, by applying new methods.

As an example of a comparatively recent development, machine learning has
begun to play a role within astronomy. Machine learning comprises a certain
subset of algorithms that can be used to find pattern in data, often involving
a training phase during which the software learns about classifying certain
kinds of data before moving on to new classifications.121212Some information
can be found in Ntampaka et al. 2019 [https://arxiv.org/abs/1902.10159].

After this general overview, let us consider astronomical data sets, and
explore ways of viewing or analyzing them. We start with astronomical images,
and a simple application software for viewing them.

##  3 SAOImage DS9 and astronomical images

When astronomy is in the public eye, a large portion of the attention goes to
spectacular astronomical images. In this section, we will take a closer look
at the scientific versions of such images. The application software we will
use in this section is called SAOImage DS9. It was developed at the
Smithsonian Astrophysical Observatory (SAO) and is available for download from
[http://ds9.si.edu/]. When you start DS9, it will look as in Fig. 33.

Figure 33: The initial screen of DS9 right after the software has started

Top left, there are fields with information about the image(s); top right, two
small windows we will talk about later. Below, there is a horizontal menu with
two rows. I will call the top row the “main menu.” Whenever you click on a
field in the top row, the bottom row will display a set of associated
commands. I will call this bottom row the “secondary menu”. Below the
secondary menu is the main image window, with a color bar (currently
greyscale) below that. There is also a horizontal menu on top — on a Mac at
the top of the screen, on Windows or on Linux at the top of the window. I will
call this the “top menu.” This is mostly a duplicate of the menu with the two
rows within the window, but it does provide some additional options, and will
occasionally be needed. Should operating systems change in the future, the
menu placement may change as well, but currently, the positioning is as I have
described. If DS9 opens additional windows, those can come with their own top
menu, different from that of the main window. Details of what the menu items
are called might vary slightly from version to version; I have used version
7.5.

###  3.1 Loading a Hubble image

Let’s load an astronomical image file. We get our file from the Hubble legacy
archive at the Space Telescope Science Institute (STScI) in Baltimore, which
operates the Hubble Space Telescopes and other space telescopes. The legacy
archive is where all the older Hubble images are stored. It can be found at
the URL [https://hla.stsci.edu/hlaview.html]. There is a helpful search field.
Let’s search for M 16, the Eagle Nebula, by entering “M16” in the search field
and pressing “Search”.

You should get a very wide result screen, the leftmost bit of which is shown
in Fig. 34.

Figure 34: Search results from the Hubble Legacy Archive

You will probably need to scroll right to see the columns 13 and 14 we are
interested in. In column 14, called “Dataset”, look for the data set names

  * •

hst_05773_05_wfpc2_f502n_wf

  * •

hst_05773_05_wfpc2_f656n_wf

  * •

hst_05773_05_wfpc2_f673n_wf

At the time of this writing, it is fairly easy to find these images. If you
don’t, you can try something different: In the search field on
[https://hla.stsci.edu/hlaview.html], click on “Advanced search”. In the main
field, enter “M16”, and in the proposal ID field, enter “05773”. The result
will be a much shorter list, including the images listed above. Even in the
far future, when the interface might have changed, searching for the proposal
ID in addition to the object name should return a list that includes those
images — even future archives should be “legacy-proof.”

The images in question were all taken on April 1, 1995, and belong to one of
the most iconic Hubble images: the pillars of creation. Each file should be
around 53 MB in size. You can download the files by either clicking the little
shopping cart icon for each image and then going to the shopping cart tab, or
by right clicking on the shopping cart icon and choosing “save link as”.

If you know astronomical abbreviations, you will be able to make some initial
sense out of these dataset names: hst, for instance is bound to mean that we
are downloading data from the Hubble Space Telescope. WFPC2 is the “Wide-Field
and Planetary Camera 2” on that telescope, wf says that we are downloading the
wide-field camera images. f502, f656 and f673 denote different filters which
have been placed before the camera for these respective images. We will
combine the three images into a colour image — but it is going to be a false-
color images, since those three filters do not correspond to red, green, and
blue!

Last but not least, when you have downloaded the images, you will notice that
the filename extension indicates that these are FITS files, in the most common
format used for scientific images in astronomy; the filename extension is
either “fits” or possibly if you are on an older Windows machine, “fit”.

###  3.2 A first look at the Eagle Nebula M16

Now that we have the files safely stored away, we can open them using DS9. To
this end, go to the main menu row and click on “file” (which is probably
highlighted to begin with); from the secondary menu that appears directly
below, choose “open”. In the usual pop-up choose-a-file window, I’ll choose
the first of the Hubble files we downloaded,
hst_05773_05_wfpc2_f502n_wf_drz.fits. Once the file is open, the DS9 window
should look as in Fig. 35.

Figure 35: Opening an HST Eagle Nebula image in DS9

The image is disappointingly black. We need to find a better brightness scale
to see what is going on. Astronomical images typically capture an amazing
dynamic range, that is, an amazing range of different brightness values for
each pixel (concretely, 65536 different brightness values per pixel, compared
with the 256 of a typical RGB pixel). Displaying such an image on a computer
monitor, or printing it, can never do this range full justice. Instead, we
need to pick and choose — which part of the brightness range do we want to
display, and which way of compressing the brightness scale shows us the most
information about the image?

There is no single right way of doing this, and there is no standard way that
will guarantee the best results for all possible astronomical pictures.
Instead, this is a matter of combining experience (your own and that of
others!) with some experimentation to arrive at a result that works for you.
You should always be aware that such a result is not a naked view of the
astronomical data — what you see is determined both by the astronomical data
and by the choices you have made in displaying that data. (Also note that
these display options do not change the image itself; the image file itself is
unchanged, and you are only changing the way you are viewing the data.)

To experiment a bit, go to the main menu and choose “scale”. From the
secondary menu that will come up, choose “zscale”. Where previously (in
“minmax” mode) the image had been displayed with the minimum pixel brightness
set to black, and the maximum pixel brightness to white, the colors are now
mapped to values closer to the median pixel brightness. As a result, your
image should look something like this in Fig. 36.

Figure 36: HST Eagle nebula image in DS9 with zscale applied

This still looks fairly raw, and rather different from the pretty astronomy
pictures you see in the media. But at least we can discern some structures.
The main window below the two rows of menus only shows part of the image. In
the overview window (second to right, on top) you can see the whole of the
image. The cyan frame in the overview window marks the part of the image that
is visible in the main window. Drag it around (left-click the mouse and drag)
to explore other parts of the image. Alternatively, you can go to “zoom” in
the main menu and choose one of the options in the secondary menu to see the
image as a whole (“zoom fit”), or in more detail (e.g. “zoom 2”).

In the “scale” menu, you can also choose another compression method instead of
linear — see how it affects the view! (Again, the image itself is not changed
by your choice.) Also, instead of the grayscale display, you can go to the
main menu point “color” and select another color map.

###  3.3 Coordinates: Navigating the image

When your cursor is on the main image, what you will see will be something
like in Fig. 37.

Figure 37: HST Eagle nebula image with cursor on a star

In this image, the tip of the cursor is placed on a star. The detailed image
(inset image top right) shows the star’s little disk, with DS9’s own
representation of the cursor on top. Let us take a closer look at the
information on the top left. “File” is simply the file name. If the header
identifies the object by name in a suitable way, that is what would be
displayed in “Object.” The “Value” field gives you the value of that
particular pixel.

Below, we have the sky position information for the cursor, given in the
“World Coordinate System” (WCS). “fk5” tells you that the sky coordinate
system is defined using the reference stars of the Fifth Fundamental Catalogue
(Fundamentalkatalog 5), which was published in 1988 by Astronomisches
Recheninstitut Heidelberg (ARI, now a part of Heidelberg University). The
coordinates themselves are those of the equatorial system, which is analogous
to latitude and longitude on Earth: longitude corresponds to the right
ascension $\alpha$ (sometimes abbreviated to RA, ra, or R.A.), latitude
corresponds to the declination $\delta$ (sometimes abbreviated to Dec, dec, or
DE).

Both RA and Dec are given here in the standard sexagesimal notation. Right
ascension is given in hours (in hour example: 18), minutes (18) and seconds
(50.2864), written as 18:18:50.2864. (In other contexts, it might be written
as $18^{h}\>18^{m}\>50^{s}.2864$, but the meaning is always the same: Each
object lies on a particular meridian (that is, on part of a great circle
through the two poles of the coordinate system). The value of the right
ascension indicates where that meridian intersects the celestial equator,
measured eastwards along the celestial equator, starting from the vernal point
— the point where the Sun crosses the celestial equator northwards, at the
time of the vernal equinox (at or around March 20 each year). For this
measurement, 24 hours correspond to the full circle. The integer hour value is
followed by the number of full minutes, with 60 minutes in a full hour, as
expected; the third position gives the seconds, again with 60 seconds in a
minute.

Declination is given in degrees northwards (positive sign) or southwards
(negative sign) from the celestial equator, which itself is at
$\delta=0^{\circ}$. In the image, the notation is again sexagesimal, with a
southward (because of the minus sign) $13$ degrees, 50 arc minutes and 2.975
arc seconds written as 13:50:2.975. (In another context, you might see this
written as $13^{\circ}50^{\prime}2^{\prime\prime}.975.$)

If you go to “WCS” in the top menu, you can change the notation from
sexagesimal to degrees; the latter will show both right ascension and
declination as a decimal number denoting degrees, in our example $274.7095$
for the right ascension, and $-13.9341$ for the declination (in both cases
with a few more significant digits).

Just like in the usual geographical coordinate system, a difference of one
degree in declination corresponds to the same length, wherever we are on the
celestial sphere, just like a difference of one degree in latitude does.
Differences in right ascension, on the other hand, correspond to smaller
angular distances the closer you go to the celestial poles, analogous to what
happens with geographic longitude. If you want to move the cursor around a bit
to estimate the angular scale of your image (how many pixels corresponding to,
say, 10 arc seconds), use declination for the purpose, not right ascension.

The ”Image” coordinates below denote the X and Y coordinate of a pixel within
the given image. In a FITS file, the pixel in the bottom left corner has the
coordinate $(1,1)$. If you zoom in, you will see that DS9 assigns this
coordinate value to the center of the pixel. The lower left corner of the
lower left corner has the pixel coordinates $(0.5,0.5)$. In our example the
“Physical” coordinate is the same as the image coordinate. If the image you
are looking at is only a part of a larger image, you are likely to find the
physical and the image coordinates differ: The physical coordinates would
still be those of the original image, $(1,1)$ the coordinates of the bottom
left pixel of the CDD camera. The image coordinates would be those of the
smaller image, $(1,1)$ the coordinates of its bottom left pixel.

How does DS9 know which sky coordinates to map to the pixel coordinates? That
(meta-)information is contained in a special part of any FITS file: the
header.

###  3.4 Meta-Information: the FITS header

Astronomical images can only be interpreted correctly if you know the
conditions under which they were produced. What filter was used? How long was
the image exposed? When and where was it taken? FITS files include this meta-
information in a dedicated section called the header. We can inspect the FITS
header of our image by going on “file” in the main menu, and choosing
“header”.

Once you do this, you are given two choices; the second one is the file name
with [SCI] appended. Start by looking at the first header, which looks like in
Fig. 38.

Figure 38: First header of the HST Eagle nebula image

You don’t need to understand this in detail (I certainly don’t!), but some
bits are fairly clear: among other things, the header lists the file name
convention, telescope and instrument name, date when this was processed, and
the filter F502N that was used to take the image. Information is encoded in a
two-part way: each particular chunk of information has a keyword (here shown
in blue, on the left) and an assigned value. The two are linked by the equal
sign. Often, this is followed by a slash, after which there is a comment with
a description of the keyword’s meaning.

If you scroll down, you can find different kinds of information: the position
of the Sun at the moment of observation is encoded there, the angle between
the Moon’s position and the pointing direction of the telescope, the
observation start time and end time. Some of those keywords are specific to
the telescope, instrument, organization or project in question. Others are
more general. Under EXPTIME, you will commonly find the exposure time in
seconds, and under DATE-OBS the date of the start of the observation in
year–month–day format. TIME-OBS gives the time of the start of the
observations; the time zone is Universal Time (UT). The multiple lines marked
HISTORY typically contain information about how the image has already been
processed, which files were used as flatfield or dark frame, and which
software was used for the processing.

The [SCI] version is more image-specific, as you can see in Fig. 39.

Figure 39: Second header of the HST Eagle nebula image

It gives you lots of details about the instrument and telescope properties
during these particular observations, and important information about the
context of the data. Here, too, some keywords will be specific to the project
in question, and others more universal. NAXIS will tell you the number of axes
you are dealing with, in this case two of them. NAXIS1 and NAXIS2 will give
you the width and height of the image, in pixels. BUNIT will give you the
units of the pixel values. In this case, it is a generic “COUNTS/S”, counts
per second, in other cases it might be Jansky per beam or similar units.

The CRVAL1, CRVAL2, CRPIX1, CRPIX2 and the CD1_1, CD1_2, CD2_1 and CD2_2
contain the information that allows you to calculate, for each pixel, the
values of the two equatorial coordinates: right ascension and declination. The
CD values correspond to a matrix mapping the two coordinate systems to each
other. Often, the axes of the WCS and the image coordinate systems are
parallel to each other. (You can check this in the overview window in the
upper right corner, second from the right, where the North and East directions
are shown as N and E, and the image coordinate directions as X and Y.) In that
case, CD1_2 and CD2_1 are zero, and for the usual square pixels, CD1_1 and
CD2_2 are equal (up to a sign) and denote the image scale. In our case, CD2_2
= 2.777777777777778E-05 denotes the degrees that correspond to the width of a
single pixel; $2.\overline{7}\cdot 10^{-5}$ degrees per pixel corresponds to
$0.1$ arc seconds per pixel. CD1_1 has an extra minus sign because the X
direction and the East direction are anti-parallel in this image.

###  3.5 Making a colour image

As a nod towards pretty astronomical pictures, let us combine partial images
into an RGB color image. Recall from section 2.1 that astronomical cameras
take black-and-white images, each through a previously chosen filter (or, in
some cases, no filter at all). These different filter images can be combined
afterwards to produce an RGB color image.

In most cases, the result will be a false-color image, however, since the
astronomical filters used will not correspond to the proper red, green and
blue filters. So let’s make a color image! Since we have played with various
switches in DS9, we should reset; the easiest way is to quit and restart the
software.

Once DS9 is up and running again, create a color frame by clicking “frame” in
the main menu, and then “rgb” in the secondary menu. A little extra menu will
pop up, which can be seen in Fig. 40.

Figure 40: The RGB popup menu

With the “red” option checked, go to “file” in the main menu and “open” in the
submenu and load hst_05773_05_wfpc2_f673n_wf_drz.fits. Next, check green and
load hst_05773_05_wfpc2_f656n_wf_drz.fits. You should already see a
superposition of the red and green parts at this stage. Finally, check blue
and load hst_05773_05_wfpc2_f502n_wf_drz.fits. There is your composite color
image, but it’s looking rather dark.

With the small RGB menu window active, go to the top menu (at the top of your
screen on a Mac, and at the top of your window frame in Windows). Under the
top menu point “Lock”, choose “scale”. That way, when you change the scale,
the change will affect all the three images equally.

By setting the scale to “sqrt” (a form of compression) and the scale limits to
99% (only possible under “scale” in the top menu) I get an image that’s pretty
close to the usual appearance of the pillars of creation, cf. Fig. 41.

Figure 41: Colour image of the “pillars of creation” with DS9

The Hubble Space Telescope images you can find in the media have had
significant extra work applied to them, up to and including a person going
over the image region by region, removing remaining impurities by hand, and
cleaning up the boundaries.

###  3.6 Catalogs

In the end, we’re not here for the pretty pictures. We’re here for the
science. And since today’s science builds upon what was done before, we will
look at ways o accessing the information that is already out there about the
region whose image we are looking at. We will work with an image from data
release 9 (DR9) of the Sloan Digital Sky Survey (SDSS). The original SDSS was
the first large-scale digital survey, a pioneering project that provided
astronomers with large amounts of data that was of consistently high quality.
Were, before, astronomers had counted themselves lucky if they could do a
statistical analysis with a few hundred galaxies, the original SDSS gave them
data on more than 900.000! The DR9 is part of the third incarnation of the
survey, SDSS-III.

I have chosen an image fairly randomly by going to
https://dr9.sdss.org/fields/ and entering the RA 20.0 and the Dec 20.0 in the
“Search by Object Coordinates” and hitting “Submit”. On the results page,
click on the link “g-band FITS” to download the file shot through the SDSS g
filter. Unzip the file to obtain the unadorned FITS file named
“frame-g-007923-5-0307.fits” which is 12.4 MB in size. Open that image with
DS9. In order to look at the image, choose a linear zscale.

Figure 42: Zoomed-out version of our SDSS image, with SDSS DR9 catalog stars
circled in green

For an image with accurate positioning data (contained, as we have seen, in
the FITS header), DS9 can show the positions of the known stars. To display
them, go to the main window’s top menu. From the dropdown menu “Analysis,” go
to the “Catalogues” dropdown sub-menu, from there to “Optical” and there,
choose “SDSS Release 9,” which is the catalog associated with the image’s data
release. If you zoom out sufficiently far, your main image will look as in
Fig. 42. The image, and the surrounding area, are covered in small green
circles, and each circle marks a star (or possibly other object) that is
listed in the SDSS DR9 catalogue.

At the same time, DS9 opens up a catalog window, as shown in Fig. 43.

Figure 43: Catalog window for the SDSS DR9 stars within our window

Near the bottom of this window is a list of objects. Each object corresponds
to one of the green circles in the image. Use the scroll bar to look at the
list’s different columns, and you will see that, in addition to the object’s
position (RA and Dec), there is information about the object’s brightness: the
different apparent magnitudes, as measured using the SDSS filter set ugriz,
corresponding to specific filters centered in the near-UV (u filter, magnitude
is given in the umag column), the blue-green part of the spectrum (gmag), the
red region (rmag), the border region between red and near-infrared (imag) and
an infrared filter beyond that (zmag).

If you click on a line in that table, corresponding to the entry for a
specific object, the object’s marker circle will blink red a few times in the
main window, and centre the view on that object In that way, you can zoom in
on specific objects and have a closer look at them. Conversely, if the main
window is active and you go to “Edit” in the top menu and, in that submenu,
choose “Catalog”, then you can click on any of the little catalog marker
circles in the image, and in the catalog window, the corresponding row will be
visible and highlighted. (Choose “None” in the “Edit” dropdown menu to go back
to the normal editing mode.)

The filter field in the catalog window allows you to filter out those parts of
the catalog that do not meet your criteria. For instance, entering “$gmag
$<17$” (where gmag is the name of the column, to which you have added a dollar
sign $, which denotes in some programming languages that this is a variable)
and clicking on the “Filter” button near the bottom of the window will only
keep those objects whose g-magnitude is smaller (brighter!) than $m_{g}=17$.
You can connect several such conditions with && for a logical “and” or $\|$
for logical “or”, e.g. “$\mbox{\$gmag}<17\;\;\|\;\;\mbox{\$umag }<16$” if any
object with either $m_{g}<17$ or $m_{u}<16$ is fine.

The sort functionality of the catalog window provides another possibility for
finding your way around the catalog data. To the right of the “Sort” marker is
a menu that lets you choose any of the columns of the catalog. To the right,
by checking a box, you can choose whether to sort the catalog using values in
that particular column in ascending or descending order (increase or
decrease).

###  3.7 Photometry with regions and statistics

Let us continue with the SDSS image we had opened in the previous section. If
you haven’t done so, choose the SDSS DR9 catalog. Sort the catalog by gmag in
increasing order so you can pick out specific values for gmag. Go to the star
with gmag 19.659. It’s at around RA 20.0714 and Dec $+19.9777$, corresponding
to the pixel coordinates $X=1819,\;Y=1215$.

Figure 44: A circular region, not selected (left) and selected (right). The
much smaller green circles are the catalog indicators and not related to the
region itself

Next, we are going to define a circular region around this star. First, I zoom
in considerably, going to “Zoom 4” and clicking “Zoom in” once. In order to do
add the region, go to the top menu and, from the dropdown menu “Edit”, choose
the item “Region.” Now we are in region mode. Still in the top menu, from the
dropdown menu “Region”, go to the sub-menu “Shapes” and select “Circle”. Now,
click on the star we have chosen. A larger green circle outline appears,
marking the region we have chosen, cf. the left part of Fig. 44. Click within
that circle, and the region is selected, indicated by the four dots framing
it, as in the right part of Fig. 44. Now you can pull the circle around by
clicking, holding and dragging. Centre it on the star.

By double clicking on the circle (while still in the “Region” edit mode), you
can bring up an extra window describing the region, cf. Fig. 45. Let us call
this the “region window,” for short.

Figure 45: The extra window for our circular region

The window allows you to read off the basic properties of your region — in the
case of a circle, center and radius. You can also change the properties by
hand, giving them specific values. We will change the radius of our circle,
making it as large as possible while keeping it sufficiently small not to
include any other stars than that in the centre. Trial-and-error suggests to
me that around 11 arc seconds is a good value in this particular case — we
just about avoid the next star to the lower left of our region. To change the
radius value to 11 arc seconds, first change the unit of the radius to
“arcsec” by choosing that value from the dropdown menu that initially reads
“Degrees”. Input the value 11 for the radius. Press “Apply” so that your
choice is applied to the region. You will probably see the region circle get a
little larger in consequence.

As long as the small region window is active, you can go to the top menu, go
to the dropdown menu “Analysis”, and choose “Statistics”. Yet another window
will pop up, this one with information about the pixels in the selected
region, cf. Fig. 46.

Figure 46: Statistics window for our region

This window provides us with several interesting bits of information. It shows
us the center of the circle, and also the pixel scale, mapping image pixels to
arc second intervals on the sky. A bit lower, we are provided with the sum of
the pixel values within the region, in this case $C_{1}=15.34$, and its
measurement error, the region’s area in square arc seconds, namely
$A_{1}=376.55$, as well as the resulting surface brightness.

The bottom line shows you additional statistical descriptors for the pixel
values within the region, namely mean, median, minimum, maximum, variance,
standard deviation and root-mean-square (rms).

You can also get DS9 to display a histogram of pixel values. Activate the
region window and, in the top menu, from the dropdown menu “Analysis,” choose
“histogram”. In our case, most of the values appear to be centered around
zero, cf. Fig. 47. That appears to be the background plus noise.

Figure 47: Histogram for the pixel values in our selected region

Let us use the statistics function for regions to perform a simple aperture
photometry measurement. In most FITS files of astronomical images, a pixel
value is proportional to the light falling onto that particular portion of the
detector (although possibly there may be an overall offset added in), that is,
proportional to the flux density of the incoming light. (When comparing
different images, you need to be a bit more circumspect — for instance, if
those images have different exposure times, and the pixel values correspond to
the energy that was collected, you will need to rescale in order to make a
valid comparison between an object in one and an object in the other image.)
In the case of our SDSS, the BUNIT keyword tells us we are dealing with a
property that is proportional to a value in Jansky, and thus a flux.

Let us make a relative brightness measurement within the image, as follows.
First, a second measurement on our first selected star, this time with a
smaller circle that just about takes in the star and its immediate
neighbourhood; I choose a radius of 8 arc seconds, and get a sum of pixel
values $C_{2}=14.17$ and area $A_{2}=200.82$. The brightness density of the
annulus bounded by the two circles, in pixel values per square arc second, is

$c_{bg}=\frac{C_{1}-C_{2}}{A_{1}-A_{2}}.$ (5)

After all, $C_{1}-C_{2}$ is the sum of the pixel values in the bigger circle
minus the sum in the smaller circle; their difference is the sum of pixel
values in the annulus. Dividing this by the area of the annulus $A_{1}-A_{2}$
in square arc seconds gives us the brightness density. In our case, we have
$c_{bg}=0.006658$ per square second.

The amount of light reaching us from within the smaller, second circle is the
light we receive from the star plus the light reaching us from the background
sources within that circle. We assume that the background contribution per
pixel from the little circle is about the same as for the surrounding annulus,
which would add up to a total background contribution of $c_{bg}\cdot A_{2}$.
Subtracting this from the sum $C_{2}$, we get an estimate for the amount of
light $l$ received from the star, namely

$l=C_{2}-\frac{C_{1}-C_{2}}{A_{1}-A_{2}}\cdot A_{2}.$ (6)

In our particular case, $l_{A}=12.83$. Note that subtracting the background
contribution also subtracts any overall offset values, so $l$ should indeed be
proportional to the amount of light that was received from this star.

Let us repeat this measurement with another star. Going down the list, I
select the one with g magnitude 17.529, near the bottom of the image frame, at
$X=1123,\;Y=34$. Repeating the same measurements with a larger circle 13 arc
seconds in radius and a smaller circle with a radius of 8 arc seconds, I
obtain $C_{1}=108.36$, $A_{1}=534.70$ square arc seconds, $C_{2}=107.57$,
$A_{2}=200.82$, and applying formula (6), I obtain $l_{B}=107.09$.

The flux ratio between the two is

$\frac{l_{B}}{l_{A}}=8.35$ (7)

so in this sense, and in this filter band, the second star is a bit over 8
times brighter than the first.

We can cross-check that against the published g magnitudes. After all, given
that the collecting area and the exposure time are the same in both cases, the
ratio of our values $l_{B}$ and $l_{A}$ should be the same as the ratio of the
intensities of the two stars. Inserting this into equation (1) for the
astronomical magnitudes, we find that we should have

$m_{A}-m_{B}=-2.5\cdot\log\left(\frac{l_{A}}{l_{B}}\right)=2.30.$ (8)

The difference in the stars’ catalog g magnitudes (gmag) is 2.13. Even with
our simple measurements, we have reproduced the relative brightness of those
stars with a deviation of 0.17 mag, corresponding to an error in their
relative flux of 17%. Acceptable for a relatively crude, off-the-cuff
measurement.

Given suitable data, we could measure changes in brightness in the same way:
Compare a non-variable star (or an ensemble of several such stars) with a star
whose brightness varies over time, and document the varying star’s light
curve. Given suitably precise measurements, you can even detect a transiting
exoplanet in this way, although this will typically require considerable
accuracy — 1% in some cases, much higher accuracy in others.

###  3.8 Profiles

DS9 also gives you the tools to create brightness profiles from astronomical
images. We will take an example image from the THINGS survey by Fabian Walter
et al.,131313A description of the THINGS survey can be found in Walter et al.
2008, https://arxiv.org/abs/0810.2125 the data of which is available for
download from the THINGS project website at
[http://www.mpia.de/THINGS/Data.html]. THINGS stands for “The HI Nearby Galaxy
Survey” and provides observations of nearby galaxies in and around the HI 21
cm line, undertaken with the Very Large Array (VLA) in Socorro, New Mexico.

For the following, please download the “ro” version of the “moment 1” file for
the galaxy NGC 3198, which is named “NGC_3198_RO_MOM1_THINGS.FITS”. This is a
radial velocity map, with each pixel showing the average speed at which atomic
hydrogen gas in that region of the galaxy moves away from us or towards us.
Open the file with DS9, and from the FITS header, you will see that the units
for each pixel are “METR/SEC,” meter per second.

Let us go into region mode again, going to the top menu, there to the dropdown
menu “Edit”, and in that menu, on “Region.” We choose a particular type of
region: in the top menu, open the dropdown menu “Region” and in the sub-menu
“Shape,” choose “Projection”.

Figure 48: A profile shape laid along the major axis of the galaxy NGC 3198’s
THINGS first moment map.

Now, click on one point of the image, hold your mouse button and drag to
another point. DS9 will join the two points with a slim rectangle and, at the
same time, open a new window that shows the profile of pixel values along that
slim rectangle. Pixel values are averaged over the direction perpendicular to
the profile; in that way, you can iron out some of the local variations
(including noise) and get an accurate representation of trends on larger
scales.141414If you want the sum of the pixel values and not the average, you
can change that in the top menu associated with the region window — under
“Analysis”, you find a “Method” sub-menus where you can set a check mark at
“Average” (default) or “Sum”. If the profile window has gotten lost, select
the region and double-click the region window; in the top menu, under
“Analysis”, uncheck and re-check the “Plot 2D”, and the profile window should
re-appear.

Click on that rectangle while still in region mode, and little dots appear —
you can drag the outer dots around to change the position, orientation and
length of the projection area, and the centre dot to change the width. In our
example, I have laid the profile, 6 pixels wide along the major axis of the
apparent ellipse formed by the galaxy NGC 3198 in the sky, see Fig. 48. The
corresponding profile curve can be seen in Fig. 49.

Figure 49: Radial velocity profile of NGC 3198 based on the first moment map.

When the profile curve window is active, you can go to the item “File” in the
top menu and, from the dropdown menu, choose “Statistic”. That will give you
minimum, maximum, mean and a few other statistical descriptors about the data.

In our case, we can see that the radial velocity grows anti-symmetrically with
distance from the galaxy’s center. The mean radial velocity is at 670 km/s,
which is the velocity at which the galaxy as a whole recedes from us. The
curve shows how different parts of the galaxy move towards us or away from us
as the galaxy rotates. (In order to determine the rotation speed of the
galaxy, we would need to de-project the radial velocities, though, as the disk
of the galaxy is tilted relative to our line of sight.)

The fact that the curve flattens out at larger radii is a piece of the puzzle
that is the evidence for dark matter: Going by the galaxy’s visible matter (in
the shape of stars and gas), one would expect the rotation speed to decrease
beyond a certain distance from the centre. The fact that it doesn’t points to
a surrounding halo of additional, non-luminous matter. Comparing the rotation
curve measurement with the expected result based on the galaxy’s brightness,
it is possible to quantify how much dark matter is needed to keep the galaxy
together. Rotation curves like this, such as the ones measured by Vera Rubin
and Kent Ford in the 1970s, are an important piece of the puzzle that has led
astronomers to postulate the existence of dark matter, which does not absorb
or emit light or other forms of electromagnetic radiation.

##  4 TOPCAT and table data

Next, let us look at application software for dealing with higher-level
astronomical data, such as catalog data or, more generally, tables. The
software is called TOPCAT, which stands for “Tool for OPerations on Catalogues
And Tables.” It was developed by the astronomer Mark Taylor at the University
of Bristol, and its tagline is “does what you want with tables”. TOPCAT can be
downloaded under this link [http://www.star.bris.ac.uk/~mbt/topcat/]. It can
do many, many things with tables, including creating diagrams, histograms,
all-sky plots and the like. We will only scratch the surface of TOPCAT’s
functionality. You can find more comprehensive tutorials on the TOPCAT page,
notably here: [http://www.star.bristol.ac.uk/~mbt/topcat/#further].

TOPCAT is based on Java, so if you haven’t got a suitably recent version of
the Java Runtime Environment (JRE) on your computer, you will need to download
it from the Oracle website at
[http://www.oracle.com/technetwork/java/index.html]. (Fortunately, installing
Java is comparatively pain-free!)

###  4.1 Opening a table file

When you open TOPCAT, it looks as shown in Fig. 50. There will be an
additional row of menu keywords — on the Windows and Linux versions at the
time of this writing (Windows), this row can be found across the top of the
window, in the current Mac operating system (I’m on macOS 10.12.6) at the top
of the screen. Again, as operating systems evolve over time, the position of
this menu might change. I will call it the top menu, for short.

TOPCAT opens additional windows for different functionality, but what Fig. 50
shows is the TOPCAT base window. If you close that window, you will close
TOPCAT.

Figure 50: The TOPCAT main window

Near the top, there is a row of icons; Fig. 51 shows a few that we will need
in the following. A mouse-over will also tell you what each icon does.

|  |
---|---|---

Load Table |  View Data |  Sky Plot
---|---|---

|  |
---|---|---

TAP |  Plane Plot |  Histogram
---|---|---

Figure 51: Some of the TOPCAT icons

Let us load an example table from a file. Specifically, go to the data page of
the Galaxy Zoo project, which can be found at [https://data.galaxyzoo.org/].
Specifically, we will look at the Galaxy zoo data release 2, table 5 “Table 5
– Main sample, spectroscopic redshifts” near the middle of the page. You can
see there are several possibilities for download: a csv file, a fits file and
a VO table file. Please download the fits file.151515Should the galaxy zoo
pages have changed drastically since the time this text was written, you can
find a local version of the same table at http://www.haus-der-
astronomie.de/working-with-astro-data The meaning of the columns can be found
under this link: [https://data.galaxyzoo.org/data/gz2/zoo2MainSpecz.txt]. The
scientific paper describing the data release is Willett et al. 2013,
[https://arxiv.org/abs/1308.3496v2]. The file itself is zipped (that is,
compressed so as to make its file size smaller). TOPCAT should be smart enough
to unzip it on its own.161616If it doesn’t, then, on Windows, you will need to
use software such as 7–zip to unpack the file. On a Linux distribution,
running the gunzip command in the command line should work (“gunzip
zoo2MainSpecz.fits.gz”). On a Mac, double clicking on the file in the finder
should do the trick.

In order to open the file, click on the “Open Table” symbol in the TOPCAT base
window (cf. Fig. 51), or go to the “File” tab in the top menu and choose “Open
Table” there. In most cases, we can leave the “Format” selection in the “Open
Table” window on “(auto)” and TOPCAT will find out on its own the type of the
table file, and open the file accordingly. If you should get an error when
attempting to open the file, you can try to specify the table type explicitly
using the “Format” dropdown menu.

Open the “Filestore Browser”, select the file choose zoo2MainSpecz.fits.gz,
the Galaxy Zoo list of classified galaxies, and click on “OK” to load it.
Topcat now looks as in Fig. 52.

Figure 52: TOPCAT after having loaded our file

The name of the table file, the number of rows and the number of columns are
displayed, and some of the icons that were initially greyed out (because they
only work once a table is loaded) are now accessible. By clicking on the “View
table data” icon (cf. Fig. 51) or by going to the top menu, choosing the
“Views” drop-down menu and clicking “Table Data”, you can get an overview of
your table data, column names included, in a new window, as in Fig. 53.

Figure 53: The TOPCAT main window with table data loaded

Click the icons next to “View table data” and you will be given either the
general table data (such as the number of columns) or an overview of the
columns and their properties, respectively. The column overview includes a
description (where available), and the type of variable stored.

###  4.2 Making a sky plot

TOPCAT provides several possibilities of viewing your data. Let’s begin with
functionality that is specific to astronomy. When you are given a catalog of
astronomical objects, you might be interested in seeing where on the celestial
sphere these objects are located. Are you looking at a full-sky survey, or a
catalog that spans only one particular region? Let us ask this question for
the Galaxy Zoo list of galaxies we have loaded into TOPCAT. Where on the
celestial sphere are these galaxies located?

To this end, go to the main window (where we clicked the “view data” icon
previously) and click the sky plot icon (cf. Fig. 51). Alternatively, in the
top horizontal menu, choose Graphics $\to$ Sky Plot. Standard catalog data
will include columns giving right ascension and declination, usually indicated
by abbreviations like RA or Dec or similar, of the objects in question. TOPCAT
is usually smart enough to find those columns, and plot your objects’
positions in the sky. The result will be displayed in an extra window, as
shown in Fig. 54.

Figure 54: The TOPCAT Skyplot window

You can click and hold with your left mouse button on the sky sphere, dragging
it around. That way, you can see that the SDSS did a detailed study of a
larger patch of the sky, but also studied selected strips of sky in nearly the
opposite direction. The physical motivation for this: If you look only in one
direction, you might miss that the universe has completely different