Lcdo. Luis Hernández

Universidad Central de Venezuela

Tarea 05

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 (Ver imágenes en este link).

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.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
import math
import statistics as stats

Traemos la imagen desde su posición dentro del repositorio

Zapatoca = plt.imread('data/zapatocaImage.jpeg')

Mostramos la imagen con sus dimensiones

plt.imshow(Zapatoca)
print("Dimensiones de la imagen:")
print(Zapatoca.shape)

Dimensiones de la imagen:
(789, 1184, 3)

Usamos la función de Matplotlib para cambiar a escala de grises

plt.imshow(Zapatoca_gray, cmap=plt.get_cmap('gray'))
plt.show()
print(Zapatoca_gray.shape)

(789, 1184)

Escogemos las estrellas usando su posición dentro de la imagen

estrellas1 = Zapatoca_gray[335:345,535:545]

plt.subplot(121)
plt.imshow(Zapatoca_gray,cmap='gray',interpolation='nearest')
plt.title('Cielo de Zapatoca') 

plt.subplot(122)
plt.imshow(estrellas1,cmap='gray',interpolation='nearest')
plt.title('Estrella 1')

plt.show()

#estrellas1

Definimos la función Gaussiana en 2 dimensiones y el error asociado a la función

#Gaussiana 2D
#p[0] amplitud ,p[1] = x0 ,p[2] = y0, p[3] = desviación ,p[4] = ctte aditiva
def gauss2d(p, x, y): 
    exponente = -((x-p[1])**2 + (y-p[2])**2) / (2*p[3]**2)
    z = p[0] * np.exp(exponente) + p[4]
    return z



#Error caso 2D
def error2d(tpl, x, y, z):
    zm = gauss2d(tpl, x, y)
    errors = zm - z
    h, w = errors.shape
    errors = errors.reshape(h*w)
    return errors

Definimos la malla que vamos a graficar en función de la estrella escogida anteriormente para luego hacer el ajuste en 2 dimensiones y luego graficar la estrella con el ajuste realizado

x = np.arange(0,10,1)
y = np.arange(0,10,1)
yy,xx = np.meshgrid(x,y)

parametrosE1 = np.array([1,5,4,1,0])

best2d1, suss = leastsq(error2d, parametrosE1, args=(xx,yy, estrellas1))
print(best2d1)

[164.57679688   3.84322132   5.37363946   1.12160047 108.72874834]

zz = gauss2d(best2d1,xx,yy)

plt.imshow(zz,cmap="gray")
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fa836f5f8d0>

Se realiza el proceso con diversas estrellas

estrellas2 = Zapatoca_gray[235:255,210:230]
plt.subplot(121)
plt.imshow(Zapatoca_gray,cmap='gray',interpolation='nearest')
plt.title('Cielo de Zapatoca') 

plt.subplot(122)
plt.imshow(estrellas2,cmap='gray',interpolation='nearest')
plt.title('Estrella 2')

plt.show()

x = np.arange(0,20,1)
y = np.arange(0,20,1)
yy2, xx2 = np.meshgrid(x,y)

parametrosE2 = np.array([1,7,8,1,0])

best2d2, suss = leastsq(error2d, parametrosE2, args=(xx2,yy2, estrellas2))
print(best2d2)

[133.32456222   8.57454845   7.88117386   1.66455007 113.28901083]

zz2 = gauss2d(best2d2,xx2,yy2)
plt.imshow(zz2,cmap="gray")
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fa8343222e8>

estrellas3 = Zapatoca_gray[135:165,380:410]
plt.subplot(121)
plt.imshow(Zapatoca_gray,cmap='gray',interpolation='nearest')
plt.title('Cielo de Zapatoca') 

plt.subplot(122)
plt.imshow(estrellas3,cmap='gray',interpolation='nearest')
plt.title('Estrella 3')

plt.show()

x = np.arange(0,30,1)
y = np.arange(0,30,1)
yy3, xx3 = np.meshgrid(x,y)

parametrosE3 = np.array([1,12,12,1,0])

best2d3, suss = leastsq(error2d, parametrosE3, args=(xx3,yy3, estrellas3))
print(best2d3)

[137.60756051  12.61832069  13.20651216   2.33491046 100.65587742]

zz3 = gauss2d(best2d3,xx3,yy3)
plt.imshow(zz3,cmap="gray")
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fa82e643940>

estrellas4 = Zapatoca_gray[300:320,470:490]
plt.subplot(121)
plt.imshow(Zapatoca_gray,cmap='gray',interpolation='nearest')
plt.title('Cielo de Zapatoca') 

plt.subplot(122)
plt.imshow(estrellas4,cmap='gray',interpolation='nearest')
plt.title('Estrella 4')

plt.show()

x = np.arange(0,20,1)
y = np.arange(0,20,1)
yy4, xx4 = np.meshgrid(x,y)

parametrosE4 = np.array([1,8,5,1,0])

best2d4, suss = leastsq(error2d, parametrosE4, args=(xx4,yy4, estrellas4))
print(best2d4)

zz4 = gauss2d(best2d4,xx4,yy4)
plt.imshow(zz4,cmap="gray")
plt.colorbar()

[162.7573849    5.55722014   8.4328912   -1.40228398 101.52830036]

<matplotlib.colorbar.Colorbar at 0x7fa836f0beb8>

estrellas5 = Zapatoca_gray[180:210,1010:1040]
plt.subplot(121)
plt.imshow(Zapatoca_gray,cmap='gray',interpolation='nearest')
plt.title('Cielo de Zapatoca') 

plt.subplot(122)
plt.imshow(estrellas5,cmap='gray',interpolation='nearest')
plt.title('Estrella 5')

plt.show()

x = np.arange(0,30,1)
y = np.arange(0,30,1)
yy5, xx5 = np.meshgrid(x,y)

parametrosE5 = np.array([1,12,11,1,0])

best2d5, suss = leastsq(error2d, parametrosE5, args=(xx5,yy5, estrellas5))
print(best2d4)

zz5 = gauss2d(best2d5,xx5,yy5)
plt.imshow(zz4,cmap="gray")
plt.colorbar()

[162.7573849    5.55722014   8.4328912   -1.40228398 101.52830036]

<matplotlib.colorbar.Colorbar at 0x7fa8341cf5f8>

se define la funcion para estudiar el parámetro de anchura a media altura FWHM

def fwhm(star, p):
    FWHM = 2 * math.sqrt(np.log(2)) * p[3]
    return FWHM

Se realiza el cáculo del parámetro FWHM y se calcula su media, mediana y desviación estándar

fwhm1 = fwhm(estrellas1, best2d1)
fwhm2 = fwhm(estrellas2, best2d2)
fwhm3 = fwhm(estrellas3, best2d3)
fwhm4 = fwhm(estrellas4, best2d4)
fwhm5 = fwhm(estrellas5, best2d5)

Vfwhm = [fwhm1, fwhm2, fwhm3, fwhm4, fwhm5]
np.array(Vfwhm)

array([ 1.86758729,  2.77165768,  3.88788094, -2.33495599, 40.38848324])

plt.title('FWHM')
plt.hist(Vfwhm, bins=10)
plt.show()
plt.clf()

<Figure size 432x288 with 0 Axes>

print("La media es:", stats.mean(Vfwhm))
print("La mediana es:", stats.median(Vfwhm))
print("La desviación estándar es:", stats.stdev(Vfwhm))

La media es: 9.316130632364123
La mediana es: 2.771657678224815
La desviación estándar es: 17.52866035304895