In [1]:
# Librerías
import numpy as np
import matplotlib.pyplot as plt
#scipy
from scipy import stats,integrate
from scipy.optimize import curve_fit
import scipy
#seaborn
import seaborn as sns
# Astropy
from astropy.io import fits
from astropy.table import QTable
from astropy.table import Table, join
import astropy.units as u
from astropy.wcs import WCS
import ligo.skymap.plot
# Cosas
import webbrowser
Funciones útiles¶
In [2]:
# Gaussiana 2d
def twoD_Gaussian(xy, amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
""" Crea una Gaussiana 2D
parámetros: Coordendas, amplitud, centrado, sigmas, angulo del eje y nivel sobre
el background"""
x, y = xy
xo = float(xo)
yo = float(yo)
a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo)
+ c*((y-yo)**2)))
return g
Normalización de Histograma y la curva de densidad¶
.- Uso Kde y normalizo a mano con el intervalo y la altura de la barra o con el intervalo y la integral de la curva
In [3]:
# Librerias
#import numpy as np
#import matplotlib.pyplot as plt
#from scipy import stats,integrate
#import scipy
#creo el vector de datos
vector=np.random.normal(size=10000)*13
# Elijo una canitdad de bins
bins=20
# Histograma con el comando de Numpy
histo_S,binEdges = np.histogram(vector,bins=bins)
# Histograma con el comando de matplotlib
out = plt.hist(vector, color='darkorange', alpha=0.5,bins=bins,edgecolor='black')
#out = plt.hist(vector,bins=bins)
# Imprimo datos parciales
#print("datos", out[0])
print("extremos del intervalo",np.max(binEdges),np.min(binEdges))
# calculo el KDE (Kernel Density Estimation)
kde1 = stats.gaussian_kde(vector)
# Considero un set de puntos para evaluar
x_eval = np.linspace(-50, 50, num=200)
# Calculo el intervalo
intervalo=(np.max(binEdges)-np.min(binEdges))/bins
print("Intervalo:",(np.max(binEdges)-np.min(binEdges))/bins)
# integral sobre la curva
print("Integral kde1",integrate.quad(kde1, -100,100))
# la suma del histograma
print("Suma histo:",np.sum(histo_S) )
# y la normalización
norma_S=np.sum(histo_S)*intervalo
print("Normalización: ",norma_S)
# Hago el dibujo multiplicando la curva normalizada a 1 por la constante calculada (norma_S)
plt.plot(x_eval,kde1(x_eval)*norma_S,'r-')
plt.show()
In [4]:
# También se podría haber hecho utilizando seaborn así:
#import seaborn as sns
sns.kdeplot(vector, fill=True,legend=True,label='Curva')
plt.legend()
plt.show()
# con Histplot (seaborn)
df = np.random.randn(100)
sns.histplot(vector, kde=True, bins=20)
plt.show()
# Con displot
#df = np.random.randn(100, 4)
sns.displot(vector, kde=True, bins=20)
plt.show()
In [5]:
# En 2D sería así:
# Librerias
#import numpy as np
#import matplotlib.pyplot as plt
#from scipy.optimize import curve_fit
# Parameters
num_points = 100000 # Number of points
mean = [0, 0] # Mean of the Gaussian distribution
cov = [[2, 0], # Covariance matrix (modificar para una elipse)
[0, 2]]
# Generate random points from a 2D Gaussian distribution
points = np.random.multivariate_normal(mean, cov, num_points)
# Scatter plot of the generated points
plt.figure(figsize=(6, 6))
plt.scatter(points[:, 0], points[:, 1], alpha=0.5, marker='.')
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.title('2D Gaussian Random Points')
plt.grid(True)
plt.show()
# Compute distances from the center
distances = np.sqrt(points[:, 0]**2 + points[:, 1]**2)
# Compute histogram with area normalization
bins = np.linspace(0, np.max(distances), 30)
hist, bin_edges = np.histogram(distances, bins=bins, density=True)
# Normalize by ring area
bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])
areas = np.pi * (bin_edges[1:]**2 - bin_edges[:-1]**2)
normalized_hist = hist / areas
# Plot density curve of distances
plt.figure(figsize=(5,5))
sns.kdeplot(distances, fill=True, color='blue', alpha=0.5)
plt.xlabel('Distance from Center')
plt.ylabel('Density')
plt.title('Density Curve of Distances from Gaussian Center')
plt.grid(True)
plt.show()
# Plot normalized density curve
plt.figure(figsize=(6, 4))
plt.plot(bin_centers, normalized_hist, marker='o', linestyle='-', color='blue')
plt.xlabel('Distance from Center')
plt.ylabel('Normalized Density')
plt.title('Area-Normalized Density Curve of Distances')
plt.grid(True)
plt.show()
# Define Gaussian function
def gaussian(x, a, mu, sigma):
return a * np.exp(-((x - mu) ** 2) / (2 * sigma ** 2))
# Fit Gaussian curve to the normalized histogram
popt, _ = curve_fit(gaussian, bin_centers, normalized_hist, p0=[1, np.mean(distances), np.std(distances)])
# Print fitted parameters
a_fit, mu_fit, sigma_fit = popt
print(f"Fitted Parameters: a = {a_fit:.4f}, mu = {mu_fit:.4f}, sigma = {sigma_fit:.4f}")
# Plot normalized density curve with Gaussian fit
plt.figure(figsize=(6, 4))
plt.plot(bin_centers, normalized_hist, marker='o', linestyle='-', color='blue', label='Data')
plt.plot(bin_centers, gaussian(bin_centers, *popt), linestyle='--', color='red', label='Gaussian Fit')
plt.xlabel('Distance from Center')
plt.ylabel('Normalized Density')
plt.title('Area-Normalized Density with Gaussian Fit')
plt.legend()
plt.grid(True)
plt.show()
# Plot normalized density curve
plt.figure(figsize=(6, 4))
plt.plot(bin_centers, hist, marker='o', linestyle='-', color='blue')
plt.xlabel('Distance from Center')
plt.ylabel('Normalized Density')
plt.title('Area-Normalized Density Curve of Distances')
plt.grid(True)
plt.show()
In [6]:
# Dibujo y ajusto la curva
#import numpy as np
#import matplotlib.pyplot as plt
#import seaborn as sns
#from scipy.optimize import curve_fit
# Parameters
num_points = 100000 # Number of points
mean = [0, 0] # Mean of the Gaussian distribution
cov = [[1, 0.5], # Covariance matrix
[0.5, 1]]
# Generate random points from a 2D Gaussian distribution
points = np.random.multivariate_normal(mean, cov, num_points)
# Scatter plot of the generated points
plt.figure(figsize=(6, 6))
plt.scatter(points[:, 0], points[:, 1], alpha=0.5, marker='o')
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.title('2D Gaussian Random Points')
plt.grid(True)
plt.show()
# Compute distances from the center
distances = np.sqrt(points[:, 0]**2 + points[:, 1]**2)
# Compute histogram with area normalization
bins = np.linspace(0, np.max(distances), 30)
hist, bin_edges = np.histogram(distances, bins=bins, density=True)
# Normalize by ring area
bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])
areas = np.pi * (bin_edges[1:]**2 - bin_edges[:-1]**2)
normalized_hist = hist / areas
# Define Gaussian function
def gaussian(x, a, mu, sigma):
return a * np.exp(-((x - mu) ** 2) / (2 * sigma ** 2))
# Fit Gaussian curve to the normalized histogram
popt, _ = curve_fit(gaussian, bin_centers, normalized_hist, p0=[1, np.mean(distances), np.std(distances)])
# Print fitted parameters
a_fit, mu_fit, sigma_fit = popt
print(f"Fitted Parameters: a = {a_fit:.4f}, mu = {mu_fit:.4f}, sigma = {sigma_fit:.4f}")
# Plot normalized density curve with Gaussian fit
plt.figure(figsize=(6, 4))
plt.plot(bin_centers, normalized_hist, marker='o', linestyle='-', color='blue', label='Data')
plt.plot(bin_centers, gaussian(bin_centers, *popt), linestyle='--', color='red', label='Gaussian Fit')
plt.xlabel('Distance from Center')
plt.ylabel('Normalized Density')
plt.title('Area-Normalized Density with Gaussian Fit')
plt.legend()
plt.grid(True)
plt.show()
Lectura de un archivo Fits que contiene una tabla Astropy¶
In [7]:
# Leo Fits y agrego los datos a una tabla
#from astropy.io import fits
#from astropy.table import QTable
#from astropy.table import Table, join
hdu_list = fits.open('J_A+A_686_A42_clusters.dat.gz.fits')
hdu_list.info()
print(hdu_list[1].columns)
data = Table(hdu_list[1].data)
data.info()
Dibujos con la proyección en el cielo¶
In [8]:
# Librerias
#import numpy as np
#from matplotlib import pyplot as plt
#import astropy.units as u
#from astropy.wcs import WCS
#import ligo.skymap.plot
#Tengo que tener ra y dec en grados
ra = [40.0,30.0,180]
dec = [-69.232,-72,0]
ax = plt.axes(projection='astro zoom',
center='4h -69d 45.4m', radius='15 deg')
ax.figure.set_size_inches(5, 5)
ax.grid()
ax.plot(ra,dec,"ro",markersize=10,transform=ax.get_transform('world'))
plt.show()
Abrir una página web en el Browser desde python¶
In [9]:
# librerías
#import webbrowser
# Sólo estos dos comandos son necesarios (o quizás uno sólo!)
#url="https://earth.google.com"
#webbrowser.open(url)
Out[9]: