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 :mod:`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:
.. code-block:: ipython3
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:
.. code-block:: ipython3
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 * =>
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 ```` es un nuevo tipo de objeto llamado `` 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.
.. code-block:: ipython3
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:
.. code-block:: ipython3
In [59]: c1.ra # Objeto
Out[59]:
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()``:
.. code-block:: ipython3
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:
.. code-block:: ipython3
In [80]: c6.galactic
In [80]:
In [81]: c6.galactic.b
Out[81]
In [82]: c_fk5 = c6.transform_to('fk5')
In [83]: c_fk5
Out[83]
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:
.. code-block:: ipython3
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
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()``:
.. code-block:: ipython3
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)
.. image:: img/staralt.png
:width: 18cm
:align: center
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:
.. code-block:: ipython3
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]:
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]:
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]:
M31
NGC1245
M81
In [241]: mags = tabla['mag']
In [242]: mags
Out[242]:
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``:
.. code-block:: ipython3
In [250]: # Objeto
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')
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``:
.. code-block:: ipython
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.
.. code-block:: bash
$ 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:
.. code-block:: ipython
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()``:
.. code-block:: ipython
In [310]: images[0].cachedataset(dir="datos")
.. code-block:: ipython
Buscamos ahora imágenes del IAC80, buscando por palabras clave:
.. code-block:: ipython
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:
.. code-block:: ipython
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:
.. code-block:: ipython
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;
.. code-block:: ipython
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:
.. code-block:: ipython
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:
.. code-block:: ipython
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:
.. code-block:: ipython
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()``:
.. code-block:: ipython
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*:
.. code-block:: ipython
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``:
.. code-block:: ipython
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:
.. code-block:: ipython
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.
.. code-block:: ipython
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.
.. code-block:: ipython
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)
.. image:: img/M81_wcsgrid.png
:width: 17cm
:align: center
Tablas FITS
-----------
Las tablas FITS funcionan muy parecido a las imágenes, con varios HDU con sus propias cabeceras y datos.
.. code-block:: ipython
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:
.. code-block:: ipython
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
^^^^^^^^^^^^^^^^^
.. code-block:: ipython
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 ```` 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.
.. code-block:: ipython
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()