Test fitting on SDSS data classified by Galaxy Zoo participants

.gitignore

@@ -0,0 +1,2 @@
+modules/__pycache__/
+

modules/galaxy_zoo.py

@@ -0,0 +1,58 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Dec  3 10:40:45 2018
+
+@author: ricktjwong
+"""
+
+from astropy.io import fits
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.image as mpimg
+from PIL import Image
+
+import sersic_profile as sp
+import mask as mk
+import process as ps
+import peak_detection as pk
+
+def rgb2gray(rgb):
+    r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]
+    gray = 0.2989 * r + 0.5870 * g + 0.1140 * b
+    return gray
+
+#hdulist = fits.open("../data/galaxyzoo1.fits")
+hdulist = fits.open("../data/psField-001740-5-0022.fit")
+data = hdulist[1].data
+print(data.dtype)
+print(len(data))
+print("Uncertain galaxy:")
+print(data[0])
+print("Spiral galaxy:")
+print(data[1])
+print("Elliptical galaxy:")
+print(data[3])
+
+#img = mpimg.imread("../data/spiral2.jpeg")
+#img = rgb2gray(img)
+#
+#plt.imshow(img)
+#
+#peaks = pk.find_peaks(img)
+#
+#total_intensity = sp.circle_intensity(img, (50, 51), 24)
+#raw, r_h = sp.half_light_radius(img, (50, 51), 24, total_intensity)
+#
+#print("Half-radius:" + str(r_h))
+#
+#y = [sp.surface_intensity(img, (50, 51), i) for i in range(1, 10)]
+#print(y)
+#x = [i for i in range(1, len(y) + 1)]
+#plt.figure()
+#plt.plot(x, y, c='b')
+#
+#I_e = sp.surface_intensity(img, (50, 51), int(r_h))
+#print(I_e)
+
+

modules/main.py

@@ -12,6 +12,7 @@
 import sersic_profile as sp
 import mask as mk
 import process as ps
+import peak_detection as pk
 
 hdulist = fits.open("../data/mosaic.fits")
 data = hdulist[0].data
@@ -30,30 +31,20 @@
 
 data = data - 3421
 y = [sp.surface_intensity(data, (2459, 3215), i) for i in range(1, 13)]
+print(y)
 x = [i for i in range(1, len(y) + 1)]
 plt.figure()
-plt.plot(x, y)
+plt.plot(x, y, c='b')
 
 I_e = sp.surface_intensity(data, (2459, 3215), int(r_h))
+print(I_e)
+print(r_h)
 
 I = [sp.I(i, I_e, int(r_h), 0.60) for i in x]
-plt.plot(x, I)
+plt.plot(x, I, "--")
 
 for j in range(1,10):
     I = [sp.I(i, I_e, int(r_h), j) for i in x]
-    plt.plot(x, I)
+    plt.plot(x, I, "--")
 
-plt.figure()
-total_intensity = sp.circle_intensity(data, (518, 3539), 3)
-raw, r_h = sp.half_light_radius(data, (518, 3539), 3, total_intensity)
-
-y = [sp.surface_intensity(data, (518, 3539), i) for i in range(1, 4)]
-x = [i for i in range(1, len(y) + 1)]
-plt.figure()
-plt.plot(x, y)
-
-I_e = sp.surface_intensity(data, (518, 3539), int(r_h))
-
-for j in range(1,10):
-    I = [sp.I(i, I_e, int(r_h), j) for i in x]
-    plt.plot(x, I)
+peaks = pk.find_peaks(data[0:1000])

modules/peak_detection.py

@@ -7,25 +7,26 @@
 
 from math import sqrt
 from skimage.feature import blob_dog, blob_log, blob_doh
-from astropy.io import fits
 import matplotlib.pyplot as plt
 
-hdulist = fits.open("../data/mosaic.fits")
-img = hdulist[0].data
 
-blobs_log = blob_log(img, max_sigma=30, num_sigma=10, threshold=.01)
-# Compute radii in the 3rd column.
-blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2)
+def find_peaks(data):
+    blobs_log = blob_log(data, min_sigma=5, max_sigma=30, num_sigma=10, threshold=.01)
+    # Compute radii in the 3rd column.
+    blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2)
+    return blobs_log
 
-fig, ax = plt.subplots()
-#coords_radius = []
-#
-for blob in blobs_log:
-    ax.set_title("LOG")
-    ax.imshow(img, interpolation='nearest')
-    y, x, r = blob
-#    coords_radius.append([x, y, r])
-    c = plt.Circle((x, y), r, color='red', linewidth=2, fill=False)
-    ax.add_patch(c)
 
-plt.show()
+def plot_blobs(data, blobs):
+    fig, ax = plt.subplots()
+    #coords_radius = []
+    #
+    for blob in blobs:
+        ax.set_title("LOG")
+        ax.imshow(data, interpolation='nearest')
+        y, x, r = blob
+    #    coords_radius.append([x, y, r])
+        c = plt.Circle((x, y), r, color='red', linewidth=2, fill=False)
+        ax.add_patch(c)
+
+    plt.show()