Wolfram Computation Meets Knowledge

Visualizing Our Place in the Milky Way Galaxy with Mathematica

In today’s world, people often forget about the wonders of the night sky. Modern conveniences provided by civilization such as electricity and lighting result in light pollution that obscures our views. Pictures like the one below that I took near Champaign, Illinois show the yellow glow of city lights that reduces the contrast with the night sky and makes it difficult to see some of the more visually stunning, but lower contrast sights like the Milky Way. But you can still make out the Milky Way in my photo as a cloudy stripe that runs up from the southern horizon during summer in the Northern hemisphere, or winter if you are in the Southern hemisphere.

Milky Way

Many people may not realize that this fuzzy band is more than just an interesting cloudy feature in the night sky. When you look at the Milky Way, you are looking directly into our home galaxy. If you think of our galaxy as a fried egg, we are sitting about 2/3 of the way from the center, out in the white part. If you look toward the yolk, your view is obscured by all the intervening white between yourself and the center of the galaxy. Although the Milky Way is not made up of eggs, the analogy helps explain the basic shape of our galaxy. It has only been in the last few decades that we have finally been able to piece together our understanding of our own galactic neighborhood and our place in it. Modern methods of astronomy, including radio and infrared observations as well as precision measurements of stellar distances from space-based observatories, have helped us put a distance scale on our models. With a bit of Mathematica code and some NASA imagery, I have reconstructed a 3D picture of our place in the Milky Way galaxy.

To start with, I need a reference point. Modern astronomy has revealed that the center of our galaxy lies at or very near a radio source called Sagittarius A* ,which is approximately 25,000 light years away. I store the distance, right ascension, and declination of Sagittarius A* as follows. Right ascension and declination are sky coordinates that are like latitude and longitude on Earth, and I’ll be using them as a reference coordinate system throughout this blog post.

sagacoords = {24841.2`3., 17.761124`8., -29.00775`6.};

I converted the above into spherical coordinates using the following transformation:


The next issue is that the plane of our galaxy is rotated compared to our equatorial coordinate system, so I need to make use of a rotation matrix to get the proper orientation. The following URL provided the necessary rotation matrix:

We can assemble the necessary rotation matrix by using the following angles along with RotationMatrix.

Rotation matrix

Now that I have a center point and a rotation matrix, I can construct a starting point. The following code uses ParametricPlot3D to construct two disks. The small red disk is a reference plane using the local equatorial coordinate system. The larger blue disk represents the plane of the Milky Way galaxy. I simply combine these two disks, but apply a rotation to the galactic disk via GeometricTransformation, and then also a translation to shift the center of the galactic disk to the location of Sagittarius A*. Of course in reality our galaxy is not perfectly flat, but this method will do for the purposes of this post.

ParametricPlot3D to construct 2 disk
2 disk image

The above is fine and dandy if all I want is a very mathematical and a kind of “sketch” view of our place in the galaxy. But I want to have something that looks more realistic. I can make use of texture mapping and AstronomicalData to take my visualization to the next step.

NASA has been able to create an illustration of our galaxy based on some of the most recent research showing that we live in a barred-spiral galaxy. I can import, crop, and rotate the image to remove extraneous black padding around the image and orient it for my purposes.

img = ImageResize[   ImageTake[    ImageRotate[initimage, 90 Degree], {656, 4944}, {656, 4944}], {256,     256}]

To represent the location of the Sun, I can use AstronomicalData to obtain its coordinates and then transform those into spherical coordinates, as we did above for Sagittarius A*.


For orientation purposes, I extract the ImageData and flip it to construct the texture.

imagedata = ImageData[img]; texture = Reverse /@ (Reverse@imagedata);

Now I can apply the texture using PlotStyle along with TextureCoordinateFunction, and also display a yellow point at the location of the Sun.

Display a yellow point at the location of the Sun

As you can see, there is an annoying black background that came along with the original image. It would be nice if I could get rid of that and make it transparent. To do this, I can manipulate the image data to append an alpha channel to the RGB data.

First, I determine the maximum RGB value, using the Norm of each RGB triple, from the extracted image data. This is needed for scaling purposes.

max = Max@Map[Norm, imagedata, {2}];

Next, I append an alpha channel to each RGB triple. The value of the alpha channel is scaled such that values near black appear totally transparent and values farther from black appear less transparent.

imagedata2 = Map[Append[#, 3 (Norm[#]/max)^3] &, imagedata, {2}];

As before, I flip the texture, but I use the modified image data in this case.

texture = Reverse /@ (Reverse@imagedata2);

Once again, I apply the texture, and this time I get a more subtle appearance with a transparent background that can be put against any color I want.

Galaxy on a transparent background

So now you can see our place in our home galaxy, but I can give the picture more value by including some other objects to provide more depth. Our galaxy is quite large, approximately 100,000 light years in diameter. On this scale, most of the stars you can see with the naked eye are very close to the Sun in this graphic, so I’ll skip those. Instead, I’ll use AstronomicalData to add a number of deep sky objects that are farther from us and are more relevant for this visualization. Open star clusters and globular star clusters should make nice additions. Globular star clusters are distributed in a halo around our galaxy, while open star clusters are primarily constrained to the galactic disk. Globular star clusters are large groups of hundreds of thousands of mainly old stars that group together in roughly spherical balls. Open star clusters are usually much younger groups of stars, usually containing far fewer stars than globular star clusters, and are found in the disk of the galaxy where most of the gas that they formed from can be found.

The Select[..., NumericQ[...]] code used below is necessary to weed out indicators of missing data, since AstronomicalData doesn’t have distance values for all objects.

Code to weed out indicators of missing data

In the following visualization, globular star clusters are red and open star clusters are green. The Sun is buried near the center of the open star clusters. None of these points are to scale, they are only used as location markers. We cannot see as far within the plane of the galaxy due to obscuring dust and gas, so this generates an observation bias making it seem like open star clusters are only found nearby. In fact, it’s only the nearby ones that are easy to see. Globular star clusters, because they are often found well above the plane our galaxy, can often be seen much farther away, since less dust and gas obscure our view.

Globular star clusters
Globular star clusters image

I created a flyby by changing the ViewPoint, and this results in even more depth.

Texture mapping images onto Mathematica graphics objects can really help visualization efforts for real-world simulations. Next time you generate Mathematica output, consider using this technique for added realism.

Download this post as a Computable Document Format (CDF) file.


Join the discussion

!Please enter your comment (at least 5 characters).

!Please enter your name.

!Please enter a valid email address.


  1. Interesting article, especially as the Subaru telescope just delivered unmatched new pictures of Andromeda (http://www.subarutelescope.org).

    There is of course an example for a 3D field of galaxies in the examples for AstronomicalData:
    equatorialToCartesian[{ra_, dec_,
    dist_}] := {dist Cos[15 ra Degree] Sin[dec Degree + Pi/2],
    dist Sin[15 ra Degree] Sin[dec Degree + Pi/2],
    dist Cos[dec Degree + Pi/2]}
    galaxies = {AstronomicalData[#, “Name”],
    AstronomicalData[#, “RightAscension”],
    AstronomicalData[#, “Declination”],
    AstronomicalData[#, “DistanceLightYears”]} & /@
    Cases[galaxies, {nam_, ra_?NumberQ, dec_?NumberQ, dist_?NumberQ} :>
    Tooltip[Sphere[equatorialToCartesian[{ra, dec, dist}],
    Scaled[0.009]], nam]]]

    It is very straightforward to use similar code for the surroundings of our Solar System. Question: Is this the most efficient code to query the curated data or do you have any suggestions to optimize it, e.g., for 1000 stars to show, possible in the correct color?

  2. The choice of which stars you want to plot can vary and when you start plotting this many objects, you need to optimize where possible. The code for querying for the data is about as optimized as it can get, but if you want to plot several thousand stars and maintain their approximate color (not taking brightness into account), the following would be an extension of your code. I’ve hard coded the PlotRange to about 2500 light years. You can make out a nearly spherical distribution of stars with the more distant ones flattening out along the galactic plane and being biased towards the blue end. This is because the stars we can see that are that far away are very luminous stars and tend to be blue. Closer to us, there tend to be more low mass red stars. These fainter red stars are harder to see that far away and so gives an observation bias favoring blue stars at further distances.
    stars = {AstronomicalData[#, “Name”],
    AstronomicalData[#, “RightAscension”],
    AstronomicalData[#, “Declination”],
    AstronomicalData[#, “DistanceLightYears”],
    AstronomicalData[#, “EffectiveTemperature”]} & /@

    goodstars =
    Cases[stars, {nam_, ra_?NumberQ, dec_?NumberQ, dist_?NumberQ,

    Plotting points is more efficient than spheres and using a single Point primitive with the VertexColors option allows for more efficient rendering.

    Point[equatorialToCartesian[#[[2 ;; 4]]] & /@ goodstars,
    VertexColors -> (ColorData[“BlackBodySpectrum”][#[[5]]] & /@
    goodstars)]}, PlotRange -> 2500, SphericalRegion -> True,
    Background -> Black]

  3. Hi, nice stuff. I am trying to figure out an efficient way to access SDSS-data directly from Mathematica, e.g. the spectra-csv data, for further analysis and processing. Couldn’t find any useful links. Thanks for your hints on this.