Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
1 result
Blame History Permalink
Forked from an inaccessible project.
Tarea-Clase-5-Ejercicio-01.md 24.92 KiB

Ejercicio Nº 1: Ejercicios para practicar numpy y optimización con scipy

 

Autor: Nicolás Fernández Cinquepalmi

 

Breve descripción: En observaciones astronómicas e imágenes en general, llamamos resolución espacial a la distancia angular minima a la que pueden estar dos fuentes puntuales de luz y aun poder ser reconocidas como objetos individuales.  

En el caso de la astronomía, este efecto tiene que ver con la dispersión de la luz al atravezar la atmósfera, la cual hace que una estrella, que debería en principio aparecer como una fuente puntual (pues las estrellas están muy lejos), aparezca en cambio como una mancha. Así, si dos estrellas están demasiado cerca sus manchas se superpondrán hasta el punto en que sea imposible distinguirlas como fuentes individuales.  

Para modelar este efecto, típicamente consideramos la acción de la atmósfera como la convolución de la imagen "perfecta" (como se vería desde el espacio) con un kernel gaussiano. El ancho de esa función gaussiana 2D caracteriza las condiciones de observación, varía con las condiciones climáticas y para cada sitio de la Tierra.  

La resolución espacial normalmente se toma como el FWHM de la gaussiana caracteristica registrada durante una observación. Es decir, si dos estrellas están a una distancia aparente en el cielo menor que el FWHM del efecto atmosférico, la luz de ambas fuentes se mezclará después de la convolución hasta el punto de impedir reconocerlas de modo individual.  

Además, la atmósfera puede interactuar de maneras distintas con la luz de distintas longitudes de onda, de manera que el ancho de la gaussiana puede ser distinto para observaciones con diferentes filtros.  

El objetivo de esta tarea es medir de forma aproximada la resolución espacial en una noche de observación en Zapatoca, Santander (Colombia), a partir de una foto del cielo estrellado.

 

Paso 1: Comenzamos con importar los módulos necesarios para este trabajo.

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import image
from matplotlib import pyplot
from PIL import Image
from scipy import optimize
from matplotlib.pyplot import figure

 

Paso 2: Cargamos la imagen desde su ubicación en el repositorio.

 

imageRGB = image.imread('/home/fernandezn/ejercicios-clase-05-datos/data/zapatocaImage.jpeg')/255
print("Dimensiones de la imagen:")
print(imageRGB.shape)
print("Dtype:")
print(imageRGB.dtype)
Dimensiones de la imagen:
(789, 1184, 3)
Dtype:
float64
pyplot.imshow(imageRGB)
pyplot.show()

png

 

Paso 3: Convertimos la imagen RGB en escala de grises haciendo un promedio de los 3 canales.

 

im_gray = np.mean(imageRGB, axis = 2)
plt.imshow(im_gray, cmap='gray')
<matplotlib.image.AxesImage at 0x7f85e058a080>

png

 

Paso 4: Realizamos cortes de la imagen, localizados en una estrella particular. Luego calculamos los centros para localizarlos en la imagen.

 

im_cut1 = im_gray[540:570,640:680]

max_index = np.where(im_cut1 >= np.max(im_cut1))
x01 = np.mean(max_index[1])
y01 = np.mean(max_index[0])

im_cut2 = im_gray[300:315,470:485]

max_index = np.where(im_cut2 >= np.max(im_cut2))
x02 = np.mean(max_index[1])
y02 = np.mean(max_index[0])

im_cut3 = im_gray[310:325,615:640]

max_index = np.where(im_cut3 >= np.max(im_cut3))
x03 = np.mean(max_index[1])
y03 = np.mean(max_index[0])

im_cut4 = im_gray[295:315,470:490]

max_index = np.where(im_cut4 >= np.max(im_cut4))
x04 = np.mean(max_index[1])
y04 = np.mean(max_index[0])

fig, axs = plt.subplots(2,2, figsize=(10,10))
axs[0,0].imshow(im_cut1, cmap='gray')
axs[0,0].plot(x01,y01,c='r',marker='+',markersize = 30)

axs[1,0].imshow(im_cut2, cmap='gray')
axs[1,0].plot(x02,y02,c='r',marker='+',markersize = 30)

axs[0,1].imshow(im_cut3, cmap='gray')
axs[0,1].plot(x03,y03,c='r',marker='+',markersize = 30)

axs[1,1].imshow(im_cut4, cmap='gray')
axs[1,1].plot(x04,y04,c='r',marker='+',markersize = 30)
[<matplotlib.lines.Line2D at 0x7f85e04be710>]

png

 

Paso 5: Definimos la función gaussiana 2D, la función de los momentos estadísticos y la función para realizar el ajuste.

 

def gaussian(height, center_x, center_y, width_x, width_y):
    """Returns a gaussian function with the given parameters"""
    width_x = float(width_x)
    width_y = float(width_y)
    return lambda x,y: height*np.exp(
                -(((center_x-x)/width_x)**2+((center_y-y)/width_y)**2)/2)

def moments(data):
    """Returns (height, x, y, width_x, width_y)
    the gaussian parameters of a 2D distribution by calculating its
    moments """
    total = data.sum()
    X, Y = np.indices(data.shape)
    x = (X*data).sum()/total
    y = (Y*data).sum()/total
    col = data[:, int(y)]
    width_x = np.sqrt(np.abs((np.arange(col.size)-x)**2*col).sum()/col.sum())
    row = data[int(x), :]
    width_y = np.sqrt(np.abs((np.arange(row.size)-y)**2*row).sum()/row.sum())
    height = data.max()
    return height, x, y, width_x, width_y

def fitgaussian(data):
    """Returns (height, x, y, width_x, width_y)
    the gaussian parameters of a 2D distribution found by a fit"""
    sigma = np.sqrt(data)
    params = moments(data)
    errorfunction = lambda p: np.ravel((gaussian(*p)(*np.indices(data.shape)) - data)/sigma)
    p, success = optimize.leastsq(errorfunction, params)
    return p

 

Paso 6: Realizamos el ajuste gaussiano para cada corte y calculamos su valor de FWHM.

 

# FIT IMAGEN 1

Xin, Yin = np.mgrid[0:30, 0:40]
data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)

plt.matshow(im_cut1, cmap=plt.cm.cubehelix)

params = fitgaussian(im_cut1)
fit = gaussian(*params)

plt.contour(fit(*np.indices(im_cut1.shape)), cmap=plt.cm.copper)
ax = plt.gca()
(height, x, y, width_x, width_y) = params

plt.text(0.95, 0.05, """
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
        fontsize=16, horizontalalignment='right',
        verticalalignment='bottom', transform=ax.transAxes,c='white')

FWHM1 = 2*np.sqrt(2*np.log(2))*width_x
print('El valor de FWHM es: ',np.round(FWHM1,2))
El valor de FWHM es:  42.97

png

# FIT IMAGEN 2

Xin, Yin = np.mgrid[0:15, 0:25]
data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)

plt.matshow(im_cut2, cmap=plt.cm.cubehelix)

params = fitgaussian(im_cut2)
fit = gaussian(*params)

plt.contour(fit(*np.indices(im_cut2.shape)), cmap=plt.cm.copper)
ax = plt.gca()
(height, x, y, width_x, width_y) = params

plt.text(0.95, 0.05, """
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
        fontsize=16, horizontalalignment='right',
        verticalalignment='bottom', transform=ax.transAxes,c='white')

FWHM2 = 2*np.sqrt(2*np.log(2))*width_x
print('El valor de FWHM es: ',np.round(FWHM2,2))
El valor de FWHM es:  23.25

png

# FIT IMAGEN 3

Xin, Yin = np.mgrid[0:14, 0:14]
data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)

plt.matshow(im_cut3, cmap=plt.cm.cubehelix)

params = fitgaussian(im_cut3)
fit = gaussian(*params)

plt.contour(fit(*np.indices(im_cut3.shape)), cmap=plt.cm.copper)
ax = plt.gca()
(height, x, y, width_x, width_y) = params

plt.text(0.95, 0.05, """
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
        fontsize=16, horizontalalignment='right',
        verticalalignment='bottom', transform=ax.transAxes,c='white')

FWHM3 = 2*np.sqrt(2*np.log(2))*width_x
print('El valor de FWHM es: ',np.round(FWHM3,2))
El valor de FWHM es:  21.82

png

# FIT IMAGEN 4

Xin, Yin = np.mgrid[0:20, 0:20]
data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)

plt.matshow(im_cut4, cmap=plt.cm.cubehelix)

params = fitgaussian(im_cut4)
fit = gaussian(*params)

plt.contour(fit(*np.indices(im_cut4.shape)), cmap=plt.cm.copper)
ax = plt.gca()
(height, x, y, width_x, width_y) = params

plt.text(0.95, 0.05, """
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
        fontsize=16, horizontalalignment='right',
        verticalalignment='bottom', transform=ax.transAxes,c='white')

FWHM4 = 2*np.sqrt(2*np.log(2))*width_x
print('El valor de FWHM es: ',np.round(FWHM4,2))
El valor de FWHM es:  37.46

png

 

Paso 7: Calculamos el valor medio, mediana y desvío estandar de los FWHM obtenidos para los cortes realizados. Por último, realizamos un histograma.

 

FWHM = [FWHM1, FWHM2, FWHM3, FWHM4]
Media = np.round(np.mean(FWHM),2)
Mediana = np.round(np.median(FWHM),2)
Desvio = np.round(np.std(FWHM),2)
print('El valor medio de FWHM es: ', Media)
print('La mediana de FWHM es: ', Mediana)
print('El desvío estándar de FWHM es: ', Desvio)

plt.hist(FWHM, bins = [0,10,20,30,40,50])
plt.title("Histograma Gray") 
plt.show()
El valor medio de FWHM es:  31.37
La mediana de FWHM es:  30.35
El desvío estándar de FWHM es:  9.07

png

 

Comentario: Cuanto menor sea el FWHM de tu fotografía mejor, ya que esto indica que la luz está más concentrada en el punto de origen que la causa y menos desparramada por los píxeles adyacentes. En este caso el valor medio es de 31.37.

 

Paso 8: Realizamos el ajuste gaussiano, de la misma forma pero ahora con el canal rojo.

 

image_red = np.copy(imageRGB) # creo una copia de la imagen para preservar la original
image_red[:,:,1]=0
image_red[:,:,2]=0
plt.title("Canal rojo")
plt.imshow(image_red)
<matplotlib.image.AxesImage at 0x7f85dcae0940>

png

im_red = np.mean(image_red, axis = 2)
red_cut1 = im_red[540:570,640:680]

max_index = np.where(red_cut1 >= np.max(red_cut1))
r_x01 = np.mean(max_index[1])
r_y01 = np.mean(max_index[0])

red_cut2 = im_red[300:315,470:485]

max_index = np.where(red_cut2 >= np.max(red_cut2))
r_x02 = np.mean(max_index[1])
r_y02 = np.mean(max_index[0])

red_cut3 = im_red[310:325,615:640]

max_index = np.where(red_cut3 >= np.max(red_cut3))
r_x03 = np.mean(max_index[1])
r_y03 = np.mean(max_index[0])

red_cut4 = im_red[295:315,470:490]

max_index = np.where(red_cut4 >= np.max(red_cut4))
r_x04 = np.mean(max_index[1])
r_y04 = np.mean(max_index[0])

fig, axs = plt.subplots(2,2, figsize=(10,10))
axs[0,0].imshow(red_cut1)
axs[0,0].plot(r_x01,r_y01,c='b',marker='+',markersize = 30)

axs[1,0].imshow(red_cut2)
axs[1,0].plot(r_x02,r_y02,c='b',marker='+',markersize = 30)

axs[0,1].imshow(red_cut3)
axs[0,1].plot(r_x03,r_y03,c='b',marker='+',markersize = 30)

axs[1,1].imshow(red_cut4)
axs[1,1].plot(r_x04,r_y04,c='b',marker='+',markersize = 30)
[<matplotlib.lines.Line2D at 0x7f85dc99b128>]

png

# FIT IMAGEN_RED 1

Xin, Yin = np.mgrid[0:30, 0:40]
data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)

plt.matshow(red_cut1, cmap=plt.cm.cubehelix)

params = fitgaussian(red_cut1)
fit = gaussian(*params)

plt.contour(fit(*np.indices(red_cut1.shape)), cmap=plt.cm.copper)
ax = plt.gca()
(height, x, y, width_x, width_y) = params

plt.text(0.95, 0.05, """
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
        fontsize=16, horizontalalignment='right',
        verticalalignment='bottom', transform=ax.transAxes,c='white')

FWHM_r1 = 2*np.sqrt(2*np.log(2))*width_x
print('El valor de FWHM es: ',np.round(FWHM_r1,2))
El valor de FWHM es:  55.27

png

# FIT IMAGEN_RED 2

Xin, Yin = np.mgrid[0:15, 0:25]
data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)

plt.matshow(red_cut2, cmap=plt.cm.cubehelix)

params = fitgaussian(red_cut2)
fit = gaussian(*params)

plt.contour(fit(*np.indices(red_cut2.shape)), cmap=plt.cm.copper)
ax = plt.gca()
(height, x, y, width_x, width_y) = params

plt.text(0.95, 0.05, """
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
        fontsize=16, horizontalalignment='right',
        verticalalignment='bottom', transform=ax.transAxes,c='white')

FWHM_r2 = 2*np.sqrt(2*np.log(2))*width_x
print('El valor de FWHM es: ',np.round(FWHM_r2,2))
El valor de FWHM es:  27.58

png

# FIT IMAGEN_RED 3

Xin, Yin = np.mgrid[0:14, 0:14]
data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)

plt.matshow(red_cut3, cmap=plt.cm.cubehelix)

params = fitgaussian(red_cut3)
fit = gaussian(*params)

plt.contour(fit(*np.indices(red_cut3.shape)), cmap=plt.cm.copper)
ax = plt.gca()
(height, x, y, width_x, width_y) = params

plt.text(0.95, 0.05, """
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
        fontsize=16, horizontalalignment='right',
        verticalalignment='bottom', transform=ax.transAxes,c='white')

FWHM_r3 = 2*np.sqrt(2*np.log(2))*width_x
print('El valor de FWHM es: ',np.round(FWHM_r3,2))
El valor de FWHM es:  26.51

png

# FIT IMAGEN_RED 4

Xin, Yin = np.mgrid[0:20, 0:20]
data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)

plt.matshow(red_cut4, cmap=plt.cm.cubehelix)

params = fitgaussian(red_cut4)
fit = gaussian(*params)

plt.contour(fit(*np.indices(red_cut4.shape)), cmap=plt.cm.copper)
ax = plt.gca()
(height, x, y, width_x, width_y) = params

plt.text(0.95, 0.05, """
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
        fontsize=16, horizontalalignment='right',
        verticalalignment='bottom', transform=ax.transAxes,c='white')

FWHM_r4 = 2*np.sqrt(2*np.log(2))*width_x
print('El valor de FWHM es: ',np.round(FWHM_r4,2))
El valor de FWHM es:  46.59

png

FWHM_r = [FWHM_r1, FWHM_r2, FWHM_r3, FWHM_r4]
Media = np.round(np.mean(FWHM_r),2)
Mediana = np.round(np.median(FWHM_r),2)
Desvio = np.round(np.std(FWHM_r),2)
print('El valor medio de FWHM es: ', Media)
print('La mediana de FWHM es: ', Mediana)
print('El desvío estándar de FWHM es: ', Desvio)

plt.hist(FWHM_r, bins = [10,20,30,40,50,60])
plt.title("Histograma Red") 
plt.show()
El valor medio de FWHM es:  38.99
La mediana de FWHM es:  37.09
El desvío estándar de FWHM es:  12.33

png

 

Comentario: El valor medio para el canal rojo es de 38.99, lo cual indica que el cielo posee altos valores de rojo.

 

Paso 9: Se realiza el mismo proceso para el canal verde.

 

image_green = np.copy(imageRGB) # creo una copia de la imagen para preservar la original
image_green[:,:,0]=0
image_green[:,:,2]=0
plt.title("Canal verde")
plt.imshow(image_green)
<matplotlib.image.AxesImage at 0x7f85dc9f7e10>

png

im_green = np.mean(image_green, axis = 2)
green_cut1 = im_green[540:570,640:680]

max_index = np.where(green_cut1 >= np.max(green_cut1))
g_x01 = np.mean(max_index[1])
g_y01 = np.mean(max_index[0])

green_cut2 = im_green[300:315,470:485]

max_index = np.where(green_cut2 >= np.max(green_cut2))
g_x02 = np.mean(max_index[1])
g_y02 = np.mean(max_index[0])

green_cut3 = im_green[310:325,615:640]

max_index = np.where(green_cut3 >= np.max(green_cut3))
g_x03 = np.mean(max_index[1])
g_y03 = np.mean(max_index[0])

green_cut4 = im_green[295:315,470:490]

max_index = np.where(green_cut4 >= np.max(green_cut4))
g_x04 = np.mean(max_index[1])
g_y04 = np.mean(max_index[0])

fig, axs = plt.subplots(2,2, figsize=(10,10))
axs[0,0].imshow(green_cut1)
axs[0,0].plot(g_x01,g_y01,c='b',marker='+',markersize = 30)

axs[1,0].imshow(green_cut2)
axs[1,0].plot(g_x02,g_y02,c='b',marker='+',markersize = 30)

axs[0,1].imshow(green_cut3)
axs[0,1].plot(g_x03,g_y03,c='b',marker='+',markersize = 30)

axs[1,1].imshow(green_cut4)
axs[1,1].plot(g_x04,g_y04,c='b',marker='+',markersize = 30)
[<matplotlib.lines.Line2D at 0x7f85dc8a5f98>]

png

# FIT IMAGEN_GREEN 1

Xin, Yin = np.mgrid[0:30, 0:40]
data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)

plt.matshow(green_cut1, cmap=plt.cm.cubehelix)

params = fitgaussian(green_cut1)
fit = gaussian(*params)

plt.contour(fit(*np.indices(green_cut1.shape)), cmap=plt.cm.copper)
ax = plt.gca()
(height, x, y, width_x, width_y) = params

plt.text(0.95, 0.05, """
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
        fontsize=16, horizontalalignment='right',
        verticalalignment='bottom', transform=ax.transAxes,c='white')

FWHM_g1 = 2*np.sqrt(2*np.log(2))*width_x
print('El valor de FWHM es: ',np.round(FWHM_g1,2))
El valor de FWHM es:  42.65

png

# FIT IMAGEN_GREEN 2

Xin, Yin = np.mgrid[0:15, 0:25]
data = gaussian(3, 100, 100, 20, 40)(Xin, Yin) + np.random.random(Xin.shape)

plt.matshow(green_cut2, cmap=plt.cm.cubehelix)

params = fitgaussian(green_cut2)
fit = gaussian(*params)

plt.contour(fit(*np.indices(green_cut2.shape)), cmap=plt.cm.copper)
ax = plt.gca()
(height, x, y, width_x, width_y) = params

plt.text(0.95, 0.05, """
x : %.1f
y : %.1f
width_x : %.1f
width_y : %.1f""" %(x, y, width_x, width_y),
        fontsize=16, horizontalalignment='right',
        verticalalignment='bottom', transform=ax.transAxes,c='white')

FWHM_g2 = 2*np.sqrt(2*np.log(2))*width_x
print('El valor de FWHM es: ',np.round(FWHM_g2,2))
El valor de FWHM es:  24.13

png