Librería astronómica astropy¶
Astropy es una librería de Python para Astronomía que incluye muchas funcionalidades comunes en astronomía. Hasta hace unos años, existían distintas librerías astronómicas independientes con funcionalidades específicas como asciidata (catálogos), APLpy (gráficos astronómicos) o pyfits (imágenes y tablas FITS), etc., que eran mantenidas por desarrolladores distintos con su propia estructura y sintaxis. En 2013 se publicó la primera versión pública de astropy que agrupa en la misma API y con clases comunes, librerías para lectura de catálogos, coordenadas, fechas y tiempo, y mucho más.
En febrero de 2018 se publicó la versión 3.0 de astropy, definiéndo la versión 2.0 como la versión con soporte a largo plazo (hasta finales de 2019) y para los usuarios de Python 2, ya que astropy v3 sólo está disponible para Python 3.
Astropy incluye librerías maduras y bien documentadas para:
Estructuras de datos y transformaciones
Constants (astropy.constants)
Units and Quantities (astropy.units)
N-dimensional datasets (astropy.nddata)
Data Tables (astropy.table)
Time and Dates (astropy.time)
Astronomical Coordinate Systems (astropy.coordinates)
World Coordinate System (astropy.wcs)
Models and Fitting (astropy.modeling)
Analytic Functions (astropy.analytic_functions)
Entrada y salida de datos
Unified file read/write interface
FITS File handling (astropy.io.fits)
ASCII Tables (astropy.io.ascii)
VOTable XML handling (astropy.io.votable)
Miscellaneous: HDF5, YAML, pickle (astropy.io.misc)
SAMP (Simple Application Messaging Protocol (astropy.samp)
Virtual Observatory Access (astropy.vo)
Otros
Cosmological Calculations (astropy.cosmology)
Convolution and filtering (astropy.convolution)
Data Visualization (astropy.visualization)
Astrostatistics Tools (astropy.stats)
Aún no hay módulos de astropy para algunas tareas como fotometría o espectroscopía, pero existen los paquetes afiliados de astropy, que aún no cumplen todas las especificaciones necesarias o no son lo suficientemente estables o maduros para ser módulos oficiales, pero son usables.
Además de ofrecer las principales utilidades astronómicas en el mismo paquete, astropy tiene la ventaja de ofrecer un framework común sobre el que crear nuevas funcionalidades usando las clases comunes que astropy ofrece.
Constantes astronómicas y unidades¶
Ya conocemos el módulo de constantes físicas de scipy, que incluye algunas constantes astronómicas pero muy pocas. astropy tiene un módulo de constantes astronómicas que lo complementa:
In [1]: from astropy import constants as C
In [2]: print(C.c)
Name = Speed of light in vacuum
Value = 299792458.0
Uncertainty = 0.0
Unit = m / s
Reference = CODATA 2014
In [3]: print(C.c.value)
299792458.0
In [4]: print(C.M_sun.value)
1.9891e+30
Este módulo se puede usar conjuntamente con el módulo units, que permite añadir unidades físicas a números y operar entre ellas:
In [5]: from astropy import units as u
In [6]: medida1 = 10.4*u.m # medida en metros
In [7]: medida2 = 245.65*u.imperial.inch # medida en pulgadas
In [8]: type(medida1) # Numero * <Unit> => <Quantity>
Out[8]: astropy.units.quantity.Quantity
In [9]: suma = medida1 + medida2 # resultado en metros
In [10]: suma = medida2 + medida1 # resultado en pulgadas
In [11]: v = (145.137*u.meter)/(3.46*u.second)
In [12]: print(v)
41.9471098266 m / s
In [13]: print(v.to('km/h'))
In [14]: print(C.c.to('pc/yr'))
0.306601393788 pc / yr
Como ve en la línea 8, la variable medida1, que es un float con una unidad <Unit> es un nuevo tipo de objeto llamado <Quantity> con solos que se pueden hacer operaciones entre sí, como suma.
Coordenadas celestes¶
El módulo coordinates ofrece clases para representar y manipular varios tipos de coordenadas. La clase principal es SkyCoord, con la que se define un objeto coordenada sobre la que se puede operar. SkyCoord admite varios formatos y opciones, que generalmente hay que indicar al definir el objeto.
In [40]: from astropy import units as u
In [41]: from astropy.coordinates import SkyCoord
In [42]: c1 = SkyCoord(10.625, 41.2, frame='icrs', unit='deg')
In [43]: c2 = SkyCoord('00h42.5m', '+41d12m')
In [44]: c3 = SkyCoord(ra=10.625*u.degree, dec=41.2*u.degree, frame='icrs')
In [45]: c4 = SkyCoord('00h42m30s', '+41d12m00s', frame='icrs')
In [46]: c5 = SkyCoord('00 42 30 +41 12 00', unit=(u.hourangle, u.deg))
In [47]: c6 = SkyCoord('00:42.5 +41:12', unit=(u.hourangle, u.deg))
In [48]: separation = c1.separation(c2) # Separacion angular de c1 respecto a c2
In [49]: print("Separation: {}".format(separation.to_string())) # separation.arcsec, separation.arcmin
Ahora que tenemos un objeto coordenadas, podemos acceder a sus valores y operar con ellos:
In [59]: c1.ra # Objeto <Longitude>
Out[59]: <Longitude 10.625 deg>
In [60]: c1.ra.radian # RA (Logitude) en radianes (float)
Out[60]: 0.18544123302439752
In [61]: c1.ra.hms # RA (Logitude) en objeto hms
Out[61]: hms_tuple(h=0.0, m=42.0, s=30.000000000000426)
Si queremos escribir coordenadas en un texto en un formato determinado, podemos usar el método to_string():
In [70]: c1.to_string('dms')
Out[70]: u'10d41m04.488s 41d16m09.012s'
In [71]: c1.to_string('decimal')
Out[71]: u'10.6846 41.2692'
In [72]: c1.to_string('hmsdms') # tranformacion a str
Out[72]: u'00h42m44.2992s +41d16m09.012s'
In [73]: c1.to_string("hmsdms", sep=":") # con separador ":"
Out[73]: u'00:42:30 +41:12:00'
La conversión de tipo de coordenadas es igual de fácil:
In [80]: c6.galactic
In [80]:
<SkyCoord (Galactic): (l, b) in deg
(121.12334339, -21.6403587)>
In [81]: c6.galactic.b
Out[81] <Latitude -21.640358696592607 deg>
In [82]: c_fk5 = c6.transform_to('fk5')
In [83]: c_fk5
Out[83]
<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
(10.62501153, 41.20000147)>
También es posible incluir distancias al crear el SkyCoord y junto con Time y EarthLocation nos permite calcular posiciones de objetos para un momento y lugar determinados. Para debemos indicar un lugar de observación con un objeto EarthLocation y la fecha en datetime o Time de astropy, además de las coordenadas del objeto:
In [90]: import numpy as np
In [91]: from astropy import units as u
In [92]: from astropy.time import Time
In [93]: from astropy.coordinates import SkyCoord, EarthLocation, AltAz, Angle
In [94]: # Resuelve el nombre por Sesame
In [95]: M44_coord = SkyCoord.from_name('M44')
In [96]: OT_lon = Angle('-16d30m35s') # Objeto <Angle>
In [97]: OT_lat = Angle('28d18m00s')
In [98]: obs_time = Time('2018-03-16 22:00:00') # Fecha de observacion
In [99]: teide_obs = EarthLocation(lat=OT_lat.degree*u.deg, lon=OT_lon.degree*u.deg,
height=2390*u.m) # Lugar de observacion en Tierra
In [100]: # El lugar de observacion también se puede obtener de la lista de observatorios de astropy (EarthLocation.get_site_names()) o dando una dirección que se resuelve con Google Maps:
In [101]: lapalma = EarthLocation.of_site("lapalma")
# Coordenadas altacimutales de M44 en el Observatorio del Teide para el día y hora ``obs_time``
In [103]: M44_AltAz = M44_coord.transform_to( AltAz(obstime=obs_time, location=teide_obs) )
In [104]: print(M44_AltAz)
In [105]: print(M44_AltAz.alt)
81d01m41.5929s
El objeto M44_AltAz contiene información de M44 para una localización y hora determinados. Podemos hacer lo mismo para un rango de fechas para conocer la posición a lo largo de la noche. Para esto primero creamos una lista o array de fechas añadiendo intervalos de tiempo con datetime.timedelta():
In [110]: import datetime
In [111]: # Intervalo de tiempo (5 minutos) para calcular las posiciones
In [112]: dt5min = datetime.timedelta(minutes=5)
In [113]: # Array con las fechas a partir de una lista de intervalos
In [114]: # obs_time es un objeto Time de astropy, pero hay que usar objetos
In [115]: # datetime para operar con fechas
In [116]: times = obs_time.datetime + dt5min*np.arange(100)
In [117]: # Array de posiciones para la lista de fechas (para la noche)
In [118]: M44_positions = M44_coord.transform_to( AltAz(obstime=times, location=teide_obs) )
In [119]: # Curva de altitud
In [120]: from matplotlib import pyplot as plt
In [120]: plt.plot(times, M44_positions.alt)
In [121]: plt.ylim(0, 90)
Tablas de datos¶
El objeto principal del módulo tables es Table, basado el arrays estructurados de numpy que permite hacer manipulaciones con datos tabulares:
In [216]: from astropy.table import Table
In [217]: names = ['M31', 'NGC1232', 'M81']
In [218]: magnitudes = [3.4, 9.9, 6.9]
In [219]: tabla = Table([names, magnitudes], names=('name', 'mag'),
meta={'name': 'Lista de galaxias'})
In [221]: tabla
Out[221]:
<Table masked=False length=3>
name mag
string56 float64
M31 3.4
NGC1245 9.9
M81 6.9
In [231]: tabla['mag'].format = '.3f' # formato de impresion por pantalla/fichero
In [232]: tabla
Out[232]:
<Table masked=False length=3>
name mag
string56 float64
M31 3.400
NGC1245 9.900
M81 6.900
In [231]: tabla['mag'].unit = 'mag' # unidades de la culumna mag
In [233]: tabla.show_in_browser(jsviewer=True) # Muestra la tabla en HTML en el navegador
In [234]: tabla.colnames
Out[234]: ['name', 'mag']
In [235]: tabla['name']
Out[235]:
<Column name='name' dtype='string56' length=3>
M31
NGC1245
M81
In [241]: mags = tabla['mag']
In [242]: mags
Out[242]:
<Column name='mag' dtype='float64' format='.3f' length=3>
3.400
9.900
6.900
In [243]: mags.data
Out[243]: array([ 3.4, 9.9, 6.9])
Se puede añadir una columna manualmente usando el método add_column, dando los valores como una lista o array o bien crear un objeto Column y añadirlo con add_columns:
In [250]: # Objeto <Column>
In [251]: distance = Column([2500000.,60990000, 11740000], name='distance', unit=u.lightyear)
In [252]: tabla.add_columns([distance])
In [253]: # Otra columna, a partir de la columna distance
In [254]: distance_kpc = Table.Column(c.to(u.kpc))
In [255]: # Hay que cambiarle el nombre, porque conserva 'distance'
In [256]: distance_kpc.name = "distance_kpc"
In [257]: tabla.add_columns([distance_pc])
In [258]: tabla['distance'] / tabla['distance_kpc'].to('lightyear')
<Quantity [ 1., 1., 1.]>
Si queremos guardar la tabla usamos write() indicando el formato incluyendo si queremos guardar en FITS, ya que los métodos de lectura y escritura de astropy están unificados y no es necesario usar astropy.fits:
In [260]: tabla.write("tabla.fits") # guarda en FITS
In [261]: tabla.write("tabla.dat", format='ascii') # guarda en ascii simple
In [262]: tabla.write("tabla.tex", format='latex') # guarda en Latex
Algunos formatos conservan las unidades en la cabecera (FITS, VOTABLE, IPAC), pero no todos.
Búsqueda en archivos astronómicos¶
pyvo es un paquete afiliado de astropy para búsquedas en el Observatorio Virtual (VO) usando los protocolos más comunes: TAP, SIAP, SSAP. Si los datos que buscamos no están en el VO, probablemente queramos usar astroquery, para consultas a catálogos generales.
$ pip3 install pyvo astroquery --user # instalar con pip
$ conda install -c pyvo astroquery # instalar con conda (Anaconda)
pyvo tiene funciones similares, pero más específicas del observatorio virtual:
In [300]: import pyvo
In [301]: from astropy.coordinates import SkyCoord
In [302]: pos = SkyCoord.from_name('M 42')
In [303]: # Busqueda de archivos de imagenes en IR
In [304]: archives = pyvo.regsearch(servicetype='image', waveband='infrared')
In [305]: for i, archive in enumerate(archives):
...: print(i, archive.res_title)
...:
In [309]: vista = archives[-1] # VISTA public surveys (VHS, VMC, VVV, VIDEO, VIKING)
In [310]: images = vista.search(pos=pos, size=0.1) # grados decimales
In [311]: images.fieldnames() # Campos disponibles
In [312]: images['Reference'] # lista de URL
Podemos descargar las imágenes con cachedataset():
In [310]: images[0].cachedataset(dir="datos")
Buscamos ahora imágenes del IAC80, buscando por palabras clave:
In [311]: # Servicio SIA del IAC80
In [312]: IAC80_sia = pyvo.regsearch(['IAC80'], servicetype='sia')
In [313]: iac80 = IAC80_sia[0] # Primera única entrada de la lista
In [314]: iac80a.search(pos=pos, size=0.1)
Astroquery permite búsquedas de datos generales y específicas para muchos archivos astronómicos, como Simbad, NED, Alma, SDSS:
In [330]: from astroquery.ned import Ned
In [331]: datos_M81 = Ned.query_object("M81")
In [331]: cerca_M81 = Ned.query_region("M81", radius=0.05 * u.deg)
En el ejemplo anterior, en lugar de buscar cerca de M81 por el nombre se puede dar un objeto SkyCoord.
Veamos ahora una búsqueda en CDS/Simbad:
In [340]: from astroquery.simbad import Simbad
In [341]: result = Simbad.query_object("GJ 569")
In [342]: print(result)
MAIN_ID RA DEC ... COO_QUAL COO_WAVELENGTH COO_BIBCODE
"h:m:s" "d:m:s" ...
———– ————- ————- … ——– ————– ——————- BD+16 2708 14 54 29.2359 +16 06 03.798 … A O 2018yCat.1345….0G
In [343]: print(result.colnames) [“MAIN_ID”, “RA”, “DEC”, “RA_PREC”, “DEC_PREC”, “COO_ERR_MAJA”, “COO_ERR_MINA”, “COO_ERR_ANGLE”, “COO_QUAL”, “COO_WAVELENGTH”, “COO_BIBCODE”]
En este caso nos devuelve campos básicos sobre el objeto, como nombre y coordenadas. Si queremos más información debemos crear una instancia específica con la información que nos interesa y también podemos quitar la que no queremos;
customSimbad = Simbad()
# ver datos disponibles
Simbad.list_votable_fields()
customSimbad.add_votable_fields('ra(ICRS)',
'dec(ICRS)',
'flux(V)',
'flux(I)',
'sptype', 'pm')
customSimbad.remove_votable_fields('coordinates')
result = customSimbad.query_object("GJ 569")
Búsqueda por patrón de nombre:
# Objetos Messier 100 a 110
messier = Simbad.query_object("M [0-1][0-9][0-9]", wildcard=True)
# Estrellas TOI
result_table = Simbad.query_object("TOI*", wildcard=True)
Trabajando con FITS¶
El subpaquete io de astropy ofrece métodos para entrada y salida de datos (lectura y escritura) de varios formatos, incluyendo tablas, XML y archivos FITS. El paquete fits es básicamente pyfits de stsci_python y permite trabajar con imágenes y tablas.
Lo primero es abrir el FITS, para luego acceder a la cabecera o los datos:
In [400]: from astropy.io import fits
In [401]: hdulist = fits.open("datos/M81_IAC80.fits")
In [402]: hdulist.info()
Out[402]:
Filename: datos/M81_IAC80.fits
No. Name Type Cards Dimensions Format
0 PRIMARY PrimaryHDU 94 (2048, 2048) float32
Con el método info() vemos la estructura del FITS, para saber cómo está contruído. Un fichero FITS está formado uno o varios HDU o Header Data Unit, cada uno de ellos puede tener su propia cabecera y datos. Es por eso que en el ejemplo anterior que la imagen de ejemplo tiene sólo un HDU (No. 0), que la primaria (PrimaryHDU), pero puede tener muchos más. Cuando consultamos una cabecera o datos, debemos entonces indicar a qué HDU nos referimos, aunque sólo haya una:
In [404]: header = hdulist[0].header
In [405]: header['OBJECT']
Out[405]: 'M81'
In [406]: header['RA']
Out[406]: '3:10:21'
Cada una de las entradas de una cabecera (cards), admite entero, string, float, etc. y se pueden modificar con la sintaxis de Python:
In [407]: header['CREATOR'] = "Ana"
In [408]: header.set('CREATOR', 'Ana')
In [409]: header['CREATOR'] = ("Ana", "Quien redujo las imagenes")
Los campos history y comment se añaden el lugar de reemplazarse.
Podemos trabajar con el objeto header como una lista (aunque no lo es exactamente). Para ver lista de nombres de la cabecera podemos usar keys():
In [410]: header[:10] # Diez primeras lineas de la cabecera
In [411]: header.keys()[:10] # Nombres de las Cards de la cabecera
Se pueden añadir campos nuevos al final como si fuesen listas con el método append o insert para hacerlo antes del campo indicado como primer parámetro, salvo que añadamos after=True, en cuyo caso lo pondrá después del campo indicado. El valor de campo se pone como segundo parámetro, que puede ser una tupla con clave, valor y comentario:
In [412]: header.insert( 'CREATOR', ('OBSERVER', 'Chio', 'Quien hizo la observacion'))
Al igual que con la cabecera, accedemos a los datos con el método data(), que devuelve un array de numpy:
In [413]: data = hdulist[0].data
In [414]: type(data)
Out[414]: numpy.ndarray
In [415]: data.shape
Out[415]: (2048, 2048)
Si queremos guardar los cambios hechos en la cabecera y datos, usamos el método writeto() para guardar en un fichero nuevo:
In [416]: hdulist.writeto("datos/M81_IAC80-new.fits", overwrite=True)
Hay que fijarse que para los cambios en la cabecera usamos la variable header y no directamente hdulist[0].header y a pesar de eso los cambios se guardan. Esto es porque header no es una copia de la cabecera, si no un puntero de memoria que apunta a hdulist[0].header; de querer realmente una copia pudimos haber hecho header2 = hdulist[0].copy().
Para sobreescribir el fichero original, debemos abrirlo en modo escritura y luego hacer flush() para pasar los cambios hechos desde el momento de apertura con open() a documento original.
In [417]: hdulist = fits.open("datos/M81_IAC80.fits", mode='update')
In [418]: hdulist[0].header['CREATOR'] = "Pepe"
In [419]: hdulist.flush() # guarda los cambios en "datos/M81_IAC80.fits
In [420]: hdulist.close() # Cierra el fichero
Con writeto() podemos indicar como segundo parámetro un array de datos, creando él una cabecera básica, o incluir un tercer parámetro con la cabecera.
In [421]: fits.writeto("datos/imagen.fits", data, header=None, overwrite=True)
World Coordinate System¶
El módulo wcs posee herramientas para trabajar con datos de WCS tomados de la cabecera de una imagen FITS y hacer transformaciones entre píxeles y coordenadas mundiales
wcs = WCS('datos/M81_IAC80.fits')
wcs.all_pix2world(1306, 575, 0)
# Pasamos de pixel a coordendas de cielo.
# El ultimo parametro se refiere al primer índice, que es 0 en numpy y C
# pero es 1 en Fortran y en ds9
coord = w.all_pix2world(1306, 575, 1)
# Coordenadas en grados
print(coords)
Ahora hacemos lo contrario, de coordenadas celestes a pixel:
pix = wcs.all_world2pix(*coord, 0)
#[array(1306.), array(575.)]
Podemos comprobar eso con las coordenadas de catálogo:
from astropy.coordinates import SkyCoord
coords = SkyCoord.from_name("M 81")
coords.to_pixel(wcs)
# (array(1052.06715444), array(999.64368801))
Es posible usar la información de astrometría en el WCS para cambiar los ejes del dibujo con imshow.
# Nuevo plot con proyeccion con la astrometria
ax = plt.subplot(projection=wcs)
plt.imshow(hdu0.data, cmap="jet", origin="bottom", vmin=300, vmax=1000)
# Malla con ejes celestes
plt.grid(color='white', ls='solid')
# Malla con ejes galacticos
gal_grid = ax.get_coords_overlay('galactic')
gal_grid.grid(color='black', ls='dotted', alpha=0.5)
Tablas FITS¶
Las tablas FITS funcionan muy parecido a las imágenes, con varios HDU con sus propias cabeceras y datos.
In [450]: hdulist = fits.open("datos/M3_photometry.fits")
In [451]: hdulist.info()
Filename: datos/M3_photometry.fits
No. Name Type Cards Dimensions Format
0 PRIMARY PrimaryHDU 14 (631,) uint8
1 BinTableHDU 21 4699R x 2C [D, D]
In [452]: tabla = hdulist[1].data
In [453]: data.shape
Out[453]: (2048, 2048)
In [454]: type(tabla)
Out[454]: function
In [455]: tabla.columns # Objetos column
In [456]: tabla.columns.names # Nombres de las columnas
In [457]: tabla.columns.info() # informacion general
In [458]: tabla.field('RAJ2000')[:5] # devuelve un array
Out[458]:
array([ 205.55550417, 205.55547917, 205.55515833, 205.55507083,
205.55503333])
Para añadir columnas nuevas o crear una tabla de cero, hay que crear cada columnas dándole la información básica:
In [470]: # Datos de nueva columna: la suma de RA y DEC (decimales)
In [471]: new_data = tabla.field('RAJ2000') + tabla.field('DEJ2000')
In [472]: # Nuevo objeto columna, con toda la informacion
In [473]: newcol = fits.Column(name='RA+DEC', format='50A', array=new_data) # 50A para string de 50 caracteres
In [474]: # Lista de definiciones de columna (ColDefs), al que anhadimos la nueva (suma de listas)
In [475]: cols = fits.ColDefs( tabla.columns + newcol)
In [477]: # Nuevo HDU
In [478]: hdu_l = fits.BinTableHDU.from_columns(cols)
In [479]: hdu_l.writeto('datos/tabla.fits')
Algunos truquitos¶
In [500]: # Lee solo la cabecera
In [501]: header = fits.getheader("datos/M81_IAC80.fits")
In [502]: # Lee solo los datos a array numpy
In [503]: data = fits.getdata("datos/M81_IAC80.fits")
In [504]: # Devuelve datos y cabecera
In [505]: data, header = fits.getdata('datos/M81_IAC80.fits', header=True)
In [506]: # Devuelve un campo de la cabecera
In [507]: filtro = fits.getval('datos/M81_IAC80.fits', 'INSFILTE')
Modelos analíticos y ajustes¶
El módulo de modelos analíticos de astropy se basa en el paquete optimize de scipy, pero ofrece modelos matemáticos para funciones comunes en Astronomía muy convenientes; además es capaz de usar las clases <Quantity> de astropy.
El módulo se basa en modelos matemáticos predefinidos.
# Gaussiana unidimensional
from astropy.modeling import models, fitting
gauss1d = models.Gaussian1D(10, 0, 7)
x = np.arange(-10, 10, 0.1)
y = gauss1d(x)
plt.plot(x, y)
O una gaussiana bidimensional:
x = np.linspace(-50, 50, 100)
y = np.linspace(-50, 50, 100)
# Arrays de parejas de coordenadas 2d de la imagen
X, Y = np.meshgrid(x, y)
z = gauss2d(X,Y)
imshow(z)
Veamos un ejemplo sencillo para hacer un ajuste gaussiano unidimensional.
In [601]: from astropy.modeling import models
# Modelo con parámetros iniciales
In [602]: g_init = models.Gaussian1D(amplitude=10, mean=0.9, stddev=0.7)
In [603]: print(g_init)
Model: Gaussian1D
Inputs: (’x’,)
Outputs: (’y’,)
Model set size: 1
Parameters:
amplitude mean stddev
--------- ---- ------
10 0.9
0.7
# creamos una distribucion normal con ruido
y = gauss1d(x) + np.random.randn(len(x))*2
# Ahora creamos una instancia de un ajustador
fitter = fitting.LevMarLSQFitter()
x = np.linspace(-50, 50, 100)
modelo = fitter(g_init, x, y)
print(ajuste)
Model: Gaussian1D
Inputs: (u’x’,)
Outputs: (u’y’,)
Model set size: 1
Parameters:
amplitude
mean
stddev
------------ ------------- --------------
10.123113704231178 -0.4248331159168961 7.122353233776502
plt.plot(x, y, 'og', label='Datos') # datos
plt.plot(x, modelo(x), 'r-', label='Modelo') # modelo
plt.legend()
La ventaja del uso de modelos es que es muy fácil crear combinaciones de modelos:
gauss1 = models.Gaussian1D(5, 0.1, 0.2)
gauss2 = models.Gaussian1D(2.0, 0.6, 0.1)
x = np.linspace(-1, 1, 200)
y = gauss1(x) + gauss2(x) + np.random.normal(0., 0.3, x.shape)
# Modelo suma de dos gaussianas
g_init = models.Gaussian1D(4, 0, 0.3) + models.Gaussian1D(2, 0.3, 0.2)
fitter = fitting.LevMarLSQFitter()
modelo = fitter(g_init, x, y)
plt.plot(x, y, '.g', label='Datos') # datos
plt.plot(x, modelo(x), 'r-', label='Modelo') # modelo
plt.legend()