Commit 853b5eac authored by rfetick's avatar rfetick
Browse files

simplify FITS writing method

parent db6ee0ce
......@@ -6,13 +6,12 @@ Created on Mon May 27 17:31:18 2019
@author: rfetick
"""
import time
import sys
import numpy as np
from scipy.fft import fft2, ifft2, fftshift, ifftshift
from astropy.io import fits
from maoppy.utils import binning
from maoppy.psfutils import moffat, gauss, oversample, moffat_center
from maoppy.psfutils import moffat, gauss, oversample, moffat_center, make_fits_hdr
__all__ = ["ParametricPSF", "ConstantPSF", "Moffat", "Gaussian",
"ParametricPSFfromPSD", "Turbulent", "Psfao"]
......@@ -48,53 +47,37 @@ class ParametricPSF:
s += "nb param: %u"%self._nparam
return s
def dict2list(dctnry):
@property
def param_name(self):
"""Ordered list of the parameters names"""
return ["param_%u"%i for i in range(self._nparam)]
def dict2list(self, dctnry):
"""Transform a dictionary of parameters into a list"""
raise NotImplementedError("This method has not been implemented for this PSF model")
return [dctnry[k] for k in self.param_name]
def list2dict(lst):
def list2dict(self, lst):
"""Transform a list of parameters into a dictionary"""
raise NotImplementedError("This method has not been implemented for this PSF model")
return {k:v for (k,v) in zip(self.param_name,lst)}
def __call__(self,*args,**kwargs):
def __call__(self, *args, **kwargs):
raise NotImplementedError("ParametricPSF is not made to be instantiated. Better use its subclasses")
def otf(self,*args,**kwargs):
def otf(self, *args, **kwargs):
"""Return the Optical Transfer Function (OTF)"""
psf = self.__call__(args,kwargs)
return fftshift(fft2(psf))
def mtf(self,*args,**kwargs):
def mtf(self, *args, **kwargs):
"""Return the Modulation Transfer Function (MTF)"""
return np.abs(self.otf(args,kwargs))
def tofits(self,param,filename,*args,keys=None,keys_comment=None,**kwargs):
def tofits(self, param, filename, *args, **kwargs):
"""Save PSF as FITS data"""
psf = self.__call__(param,*args,**kwargs)
hdr = _getfitshdr(param,keys=keys,keys_comment=keys_comment)
psf = self.__call__(param, *args, **kwargs)
hdr = make_fits_hdr(param, keys=self.param_name)
hdu = fits.PrimaryHDU(psf, hdr)
hdu.writeto(filename, overwrite=True)
def _getfitshdr(param,keys=None,keys_comment=None):
"""Define a header to save a .fits file"""
if keys is None:
keys = ["PARAM %u"%(i+1) for i in range(len(param))]
if len(keys)!=len(param):
raise ValueError("When defined, `keys` must be same size as `param`")
if keys_comment is not None:
if len(keys_comment)!=len(param):
raise ValueError("When defined, `keys_comment` must be same size as `param`")
hdr = fits.Header()
hdr["HIERARCH ORIGIN"] = "MAOPPY automatic header"
hdr["HIERARCH CREATION"] = (time.ctime(),"Date of file creation")
for i in range(len(param)):
if keys_comment is None:
hdr["HIERARCH PARAM "+keys[i]] = param[i]
else:
hdr["HIERARCH PARAM "+keys[i]] = (param[i],keys_comment[i])
return hdr
hdu.writeto(filename) # overwrite ?
class ConstantPSF(ParametricPSF):
......@@ -122,6 +105,11 @@ class Moffat(ParametricPSF):
self._xxyy[0] -= npix[0]//2
self._xxyy[1] -= npix[1]//2
@property
def param_name(self):
"""Ordered list of the parameters names"""
return ["alpha_x","alpha_y","theta","beta"]
def __call__(self, param, dx=0, dy=0):
"""
Parameters
......@@ -134,11 +122,6 @@ class Moffat(ParametricPSF):
"""
paramdxdy = np.concatenate((param,[dx,dy]))
return moffat(self._xxyy, paramdxdy, norm=self.norm)
def tofits(self,param,filename,*args,keys=None,**kwargs):
if keys is None:
keys = ["ALPHA_X","ALPHA_Y","THETA","BETA"]
super().tofits(param,filename,*args,keys=keys,**kwargs)
class Gaussian(ParametricPSF):
......@@ -153,6 +136,11 @@ class Gaussian(ParametricPSF):
self._xxyy[1] -= npix[0]/2
self._xxyy[0] -= npix[1]/2
@property
def param_name(self):
"""Ordered list of the parameters names"""
return ["sigma_x","sigma_y","theta"]
def __call__(self, param, dx=0, dy=0):
"""
Parameters
......@@ -164,11 +152,6 @@ class Gaussian(ParametricPSF):
"""
paramdxdy = np.concatenate((param,[dx,dy]))
return gauss(self._xxyy, paramdxdy)
def tofits(self,param,filename,*args,keys=None,**kwargs):
if keys is None:
keys = ["SIGMA_X","SIGMA_Y","THETA"]
super().tofits(param,filename,*args,keys=keys,**kwargs)
#%% MASTER CLASS
......@@ -501,6 +484,11 @@ class Turbulent(ParametricPSFfromPSD):
bounds_up = [np.inf, np.inf]
self.bounds = (bounds_down,bounds_up)
@property
def param_name(self):
"""Ordered list of the parameters names"""
return ["r0","lext"]
def psd(self, parampsd, grad=False):
"""Compute the PSD model from parameters
PSD is given in [rad²/f²] = [rad² m²]
......@@ -623,16 +611,6 @@ class Psfao(ParametricPSFfromPSD):
return ["r0","bck","amp","alpha","ratio","theta","beta"]
def dict2list(self, dctnry):
"""Transform a dictionary of parameters into a list"""
return [dctnry[k] for k in self.param_name]
def list2dict(self, lst):
"""Transform a list of parameters into a dictionary"""
return {k:v for (k,v) in zip(self.param_name,lst)}
def psd(self, parampsd, grad=False, fovnorm=False):
"""Compute the PSD of the electromagnetic phase
PSD is given in [rad²/f²] = [rad² m²]
......@@ -733,20 +711,18 @@ class Psfao(ParametricPSFfromPSD):
return psd, integral
def tofits(self, param, filename, *args, keys=None, overwrite=False, **kwargs):
if keys is None:
keys = self.param_name
keys_comment = ["Fried parameter [m]",
"PSD AO area constant C [rad2]",
"PSD AO area Moffat variance A [rad2]",
"PSD AO area Moffat alpha [1/m]",
"PSD AO area Moffat sqrt(ax/ay) ratio",
"PSD AO area Moffat theta [rad]",
"PSD AO area Moffat beta"]
def tofits(self, param, filename, *args, **kwargs):
keys_comment = ["Fried parameter [m]",
"AO background [rad2 m2]",
"AO Moffat variance A [rad2]",
"AO Moffat alpha [1/m]",
"AO Moffat sqrt(ax/ay) ratio",
"AO Moffat theta [rad]",
"AO Moffat beta"]
# redefine tofits() because extra hdr is required
psf = self.__call__(param,*args,**kwargs)
hdr = _getfitshdr(param,keys=keys,keys_comment=keys_comment)
psf = self.__call__(param, *args, **kwargs)
hdr = make_fits_hdr(param, keys=self.param_name, keys_comment=keys_comment)
hdr["HIERARCH SYSTEM"] = (self.system.name,"System name")
hdr["HIERARCH SYSTEM D"] = (self.system.D,"Primary mirror diameter")
......@@ -755,4 +731,4 @@ class Psfao(ParametricPSFfromPSD):
hdr["HIERARCH LEXT"] = (self.Lext,"Von-Karman outer scale")
hdu = fits.PrimaryHDU(psf, hdr)
hdu.writeto(filename, overwrite=overwrite)
hdu.writeto(filename)
......@@ -6,10 +6,10 @@ Created on Tue Jan 18 12:02:10 2022
"""
import numpy as np
from astropy.io import fits
__all__ = ["oversample", "reduced_coord", "reduced_center_coord", "moffat",
"moffat_center", "gauss"]
"moffat_center", "gauss", "make_fits_hdr"]
def oversample(samp, fixed_k = None):
......@@ -214,4 +214,22 @@ def gauss(xxyy,param):
ax = np.sqrt(2)*param[0]
ay = np.sqrt(2)*param[1]
uu = reduced_coord(xxyy,ax,ay,param[2],param[3],param[4])
return np.exp(-uu) / (2*np.pi*param[0]*param[1])
\ No newline at end of file
return np.exp(-uu) / (2*np.pi*param[0]*param[1])
def make_fits_hdr(param, keys=None, keys_comment=None):
"""Define a header to save a .fits file"""
if len(keys)!=len(param):
raise ValueError("`keys` must be same size as `param`")
if keys_comment is not None:
if len(keys_comment)!=len(param):
raise ValueError("When defined, `keys_comment` must be same size as `param`")
hdr = fits.Header()
hdr["HIERARCH ORIGIN"] = "MAOPPY automatic header"
for i in range(len(param)):
if keys_comment is None:
hdr["HIERARCH PARAM "+keys[i]] = param[i]
else:
hdr["HIERARCH PARAM "+keys[i]] = (param[i],keys_comment[i])
return hdr
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment