Getting Hands Dirty in Python

I learned python about a year and a half ago and it's been a huge time saver. I've moved from a language where you had to write a library for everything to a language where a library for everything already exists. Here's something I wrote quickly this morning. Very easy, and you can do it too. If you're interested in better understanding your astronomy data for educational purposes, this post will help get you started.
Before we begin, you'll need to have installed these libraries:
Python
pylab
rawkit
Pillow
LibRaw (Without this library, we couldn't do this thanks to developers for all this hard work!)
I'm happy to provide more instructions if need be.
I finally bought a mount (Celestron AVX) and I'm beginning to take my own images.



Here is some quick code to do this. Note, some steps are not very 'pythonic'. I'm a strong believer in breaking conventions to better understand data, and worry about conforming to them when packaging code.
This allows for much quicker analysis and gaining a deeper understanding of the data, which is *much* more important than the code!!! (coding is easy, making sense of data is the challenge...)

In this code, you'll learn to:
1. Extract the raw data from your image
2. Extract each color from your image (if needed), and make a grey scale image
3. Rescale your image in different ways to try to better emphasize features without brining out the noise

This code should be fairly straightforward and is meant just as a motivational exercise. There are better ways to do these things, and nice statistics you can use to better understand your data (like histograms). I hope this motivates you in the same way it excites me to have so much control over my data!

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99

''' Astronomy grey image analysis ''' # This imports functions into global namespace # useful for quick hacking, NOT recommended for writing programs from pylab import * # makes plots appear if in interactive mode (comment out if not) ion() # used for saving images from PIL import Image # used for reading Raw images from rawkit.raw import Raw # choose your filename here import os.path filename = os.path.expanduser("~/pictures/17-03-18-custer-whirlpool/IMG_3787.CR2") # Grab the data, really easy im = Raw(filename) data = np.array(im.raw_image()) # The data is in a Bayer Color Matrix here for Canon, so let's # just quickly extract each color by choosing the right pixels # quick way to get from Bayer matrix. Note there are way more # sophisticated algorithms out there than this method... But this # yields the rawest untouched data possible # # On Rebel T3, pixels are: # [ R G ] # [ G B ] # so index every other row/column from the right start (0 or 1) # start:end:skip means start at start, end to end and skip every other row # if "end" not specified, then just go to end of array # so 0::2 means start at 0, go to the end, and skip every other row/column # data is indexed as data[row, colum]a (row goes up/down in image, column goes # left/right) R = data[0::2,0::2] # for green, you average the two together G = (data[0::2,1::2]+data[1::2,0::2])*.5 B = data[1::2,1::2] # sum colors to get greyscale (note this isn't 100% correct, you need # to weight colors appropriately, but this is just messing around first) Grey = R + G + B # take logarithm, sometimes this could look better LGrey = np.log10(Grey) # try other combinations Grey_squared = Grey**2 # so now, let's plot the image, plot it first without "vmin" and "vmax" # and mouse over the plot to see what values you see # when happy, choose a good "vmin" and "vmax", this worked for this image set_cmap("Greys_r") figure(4);clf(); imshow(LGrey[200:900,600:1300],vmin=3.8,vmax=3.9) figure(5);clf(); imshow(Grey[200:900,600:1300],vmin=10**3.8,vmax=10**3.9) figure(6);clf(); imshow(Grey_squared[200:900,600:1300], vmin=4e7, vmax=6e7) # now, if we want to save image, it needs to have its pixel values go from # 0 to 2**8-1 (255). You don't want a larger range since your eye can't really tell # the difference beyond that, roughly # normalizing function: def normalize(img, mn, mx, logimg=True): ''' normalize to a uint8 This is also known as byte scaling. ''' dynamic_range = 2**8-1 img = img-mn img /= (mx-mn) img *= dynamic_range img = np.minimum(np.maximum(img, 0), dynamic_range) img = img.astype(np.uint8) return img img_grey = normalize(Grey[200:900, 600:1300], 10**3.8, 10**3.9) limg_grey = normalize(LGrey[200:900, 600:1300], 3.8, 3.9) img_grey_squared = normalize(LGrey[200:900, 600:1300], 4e7, 6e7) # Image comes from the PIL library (imported at top) img = Image.fromarray(img_grey) limg = Image.fromarray(limg_grey) img_grey_squared = Image.fromarray(img_grey_squared) # saving img.save("../storage/01-whirlpool-oneimage.png") limg.save("../storage/01-whirlpool-log-oneimage.png") limg.save("../storage/01-whirlpool-squared-oneimage.png")

Finally, you can run the code by running it in ipython (interactive mode):

$ ipython
[1] : %run mycode.py

or just python:

$ python run mycode.py

The whirpool galaxy. 5 min exposure taken with a Canon Rebel T3 on an Orion Maksutov-Cass 127mm on a Celestron AVX mount, autoguided.

The whirpool galaxy, logarithm of image. 5 min exposure taken with a Canon Rebel T3 on an Orion Maksutov-Cass 127mm on a Celestron AVX mount, autoguided.

The whirpool galaxy, square of image. 5 min exposure taken with a Canon Rebel T3 on an Orion Maksutov-Cass 127mm on a Celestron AVX mount, autoguided.

Right now, all images look the same, but there will be cases where different transformations help (this is similar to playing with the linearization curve). We'll mess around with colors later.