Skip to content

GitLab

Commit 45eeaff3 authored by Vigeesh Gangadharan's avatar Vigeesh Gangadharan
Browse files

works with single map

parent c3b61987
Pipeline #2910 passed with stage
in 5 minutes and 4 seconds
Hide whitespace changes
Inline Side-by-side
grisinv/invert.py
View file @ 45eeaff3
......@@ -212,6 +212,9 @@ def vfisv(
vfisv_data = VFISVpackage(path, id, line, width)
print(vfisv_data.SI.dtype,vfisv_data.line_props,vfisv_data.map_props)
#sys.exit()
# Spawn jobs to the grisinv:master+slaves
pycomm = MPI.COMM_SELF.Spawn(get_executable(), maxprocs=numproc)
......@@ -279,6 +282,8 @@ def vfisv(
data, sifit, sqfit, sufit, svfit, vfisv_data
)
print('here',np.shape(inv_data['inte']),np.shape(stokes_data['Stokes_I']))
# change the vfisv_data to a header list here and pass the header for return
header = make_header(vfisv_data)
......@@ -330,7 +335,11 @@ class FitsCardGroup:
def get_imap_list(vfisv_params):
"""Get the imap number as a list corresponding to the sorted filenames"""
return (vfisv_params.imap - 1).astype("int")
if vfisv_params.header['INSTRUME'] == 'GRIS':
return (vfisv_params.imap - 1).astype("int")
elif vfisv_params.header['INSTRUME'] == 'GRIS-IFU':
return int(vfisv_params.imap - 1)#.astype("int")
def get_istep_list(vfisv_params):
......@@ -338,7 +347,7 @@ def get_istep_list(vfisv_params):
return (vfisv_params.istep - 1).astype("int")
def get_result_shape(vfisv_params):
def get_result_shape_slit(vfisv_params):
"""Get the shape of the data"""
nmaps = vfisv_params.header["NMAPS"]
......@@ -348,6 +357,17 @@ def get_result_shape(vfisv_params):
return [nmaps, nlength, nsteps]
def get_result_shape_ifu(vfisv_params):
"""Get the shape of the data"""
nmaps = vfisv_params.header["NMAPS"]
nsteps = vfisv_params.header["NAXIS1"]
nlength = vfisv_params.header["NAXIS2"]
return [nmaps, nlength, nsteps]
def get_stokes_shape(vfisv_params):
"""Get the shape of the data"""
......@@ -383,9 +403,6 @@ def get_quantity(data, vfisv_params, index: int, units=Unit.dimensionless_unscal
the array in the proper shape and units
"""
# Get the shape for the output
data_shape = get_result_shape(vfisv_params)
# Get the "imap" ordering
imaps = get_imap_list(vfisv_params)
......@@ -395,12 +412,25 @@ def get_quantity(data, vfisv_params, index: int, units=Unit.dimensionless_unscal
# the lx_range is the range the filenames are ordered.
lx_range = np.arange(vfisv_params.LX)
# initialize the empty numpy array
quantity = np.zeros(data_shape)
# reorder according to imap and istep.
# TODO: Check for the flip that is caused while reading the fits file
quantity[imaps, :, isteps] = data[lx_range, ::-1, index]
if vfisv_params.header['INSTRUME'] == 'GRIS':
# Get the shape for the output
data_shape = get_result_shape_slit(vfisv_params)
# initialize the empty numpy array
quantity = np.zeros(data_shape)
# print('GRIS')
quantity[imaps, :, isteps] = data[lx_range, ::-1, index]
elif vfisv_params.header['INSTRUME'] == 'GRIS-IFU':
# print('GRIS-IFU')
data_shape = get_result_shape_ifu(vfisv_params)
# initialize the empty numpy array
quantity = np.zeros(data_shape)
quantity[imaps, :, :] = data[:, :, index].T
#quantity[imaps, :, isteps] = data[lx_range, ::-1, index]
#print('here',data_shape, np.shape(data),np.shape(quantity))
return quantity * units
......@@ -425,9 +455,6 @@ def get_stokes_fits(data, vfisv_params):
the array in the proper shape
"""
# Get the shape for the output
data_shape = get_stokes_shape(vfisv_params)
# Get the "imap" ordering
imaps = get_imap_list(vfisv_params)
......@@ -437,12 +464,22 @@ def get_stokes_fits(data, vfisv_params):
# the lx_range is the range the filenames are ordered.
lx_range = np.arange(vfisv_params.LX)
# initialize the empty numpy array
quantity = np.zeros(data_shape)
# reorder according to imap and istep.
# TODO: Check for the flip that is caused while reading the fits file
quantity[imaps, :, isteps, :] = data[lx_range, ::-1, :]
# quantity[imaps, :, isteps, :] = data[lx_range, ::-1, :]
if vfisv_params.header['INSTRUME'] == 'GRIS':
# Get the shape for the output
data_shape = get_stokes_shape(vfisv_params)
# initialize the empty numpy array
quantity = np.zeros(data_shape)
# print('GRIS')
quantity[imaps, :, isteps, :] = data[lx_range, ::-1, :]
elif vfisv_params.header['INSTRUME'] == 'GRIS-IFU':
# print('GRIS-IFU')
data_shape = get_result_shape_ifu(vfisv_params)
# initialize the empty numpy array
quantity = np.zeros(data_shape)
quantity = data[:, ::-1, :]
return quantity
......@@ -1564,19 +1601,39 @@ class VFISVpackage:
self.filenames = None # list containing filenames
# Populate all parameters
# Get the list of filenames
self.get_filenames()
# Get the pixels that define the boundary
self.get_pixel_bounds()
self.get_data()
#get data depending on slit or IFU
if self.header['INSTRUME']=='GRIS':
#print('GRIS')
self.get_data_slit()
elif self.header['INSTRUME']=='GRIS-IFU':
#print('GRIS-IFU')
self.get_data_ifu()
self.get_map_props()
self.get_line_props()
self.get_rel_wave()
self.nmaps = self.header["NMAPS"] # Number of maps
self.nsteps = self.header["NSTEPS"] # Number of steps per map
self.LX = len(self.filenames) # Number of slits
if self.header['INSTRUME']=='GRIS':
#print('GRIS')
self.LX = len(self.filenames) # Number of slits
elif self.header['INSTRUME']=='GRIS-IFU':
#print('GRIS-IFU')
self.LX = self.header["NAXIS1"] # width
self.LY = self.header["NAXIS2"] # Length of the slit
self.LZ = 16 # VFISV data structure
print(self.LX,self.LY)
@staticmethod
def get_line_list():
"""Get the list of available lines"""
......@@ -1627,11 +1684,18 @@ class VFISVpackage:
def get_map_props(self):
"""Create the integer parameters buffer"""
if self.header['INSTRUME']=='GRIS':
#print('GRIS')
LX = len(self.filenames)
elif self.header['INSTRUME']=='GRIS-IFU':
#print('GRIS-IFU')
LX = self.header['NAXIS1']
self.map_props = np.array(
(
1, # YMIN (for VFISV)
self.header["NAXIS2"], # YMAX
len(self.filenames), # LX
LX, # LX
self.header["NAXIS2"], # LY
self.pix_end - self.pix_ini + 1, # NUMW
1, # LINI
......@@ -1768,7 +1832,7 @@ class VFISVpackage:
if self.line == float(iline.split(",")[1]):
self.geff = iline.split(",")[3]
def get_data(self):
def get_data_slit(self):
"""
Get the stokes vector from the L1 data
"""
......@@ -1807,7 +1871,7 @@ class VFISVpackage:
self.imap[iff],
self.istep[iff],
self.time[iff],
) = self.read_fits(self.filenames[iff])
) = self.read_fits_slit(self.filenames[iff])
# Switch to a FORTRAN array
self.SI = np.asfortranarray(self.SI)
......@@ -1820,6 +1884,82 @@ class VFISVpackage:
# Truncate unfullfilled maps
self.truncate_maps()
def get_data_ifu(self):
"""
Get the stokes vector from the L1 data
"""
nx = self.header["NAXIS1"]
ny = self.header["NAXIS2"]
# Initialize the Stokes vector and continuum intensity arrays for VFISV
self.SI = np.zeros(
(len(self.filenames),nx, ny, self.pix_end - self.pix_ini + 1),
dtype="float32",
)
self.SQ = np.zeros_like(self.SI, dtype="float32")
self.SU = np.zeros_like(self.SI, dtype="float32")
self.SV = np.zeros_like(self.SI, dtype="float32")
self.ICONT = np.zeros([len(self.filenames),nx, ny], dtype="float32")
self.imap = np.ones(len(self.filenames), dtype="int")
self.istep = np.zeros(len(self.filenames), dtype="int")
self.time = np.zeros(len(self.filenames), dtype="datetime64[s]")
# # Iterate through the L1 split files/maps
# # for iff in tqdm(
# # range(len(self.filenames)),
# # bar_format=termcolors.text
# # + "Reading fits files {n_fmt}/{total_fmt}: {desc}|{bar:30}|"
# # + termcolors.end
# # + "{percentage:3.0f}% ({elapsed_s:3.1f}s)",
# # ):
#
# (
# self.SI,
# self.SQ,
# self.SU,
# self.SV,
# self.ICONT,
# self.imap,
# #self.istep[iff],
# self.time,
# ) = self.read_fits_ifu(self.filenames[0])
#
for iff in tqdm(
range(len(self.filenames)),
bar_format=termcolors.text
+ "Reading fits files {n_fmt}/{total_fmt}: {desc}|{bar:30}|"
+ termcolors.end
+ "{percentage:3.0f}% ({elapsed_s:3.1f}s)",
):
(
self.SI[iff, ...],
self.SQ[iff, ...],
self.SU[iff, ...],
self.SV[iff, ...],
self.ICONT[iff, :],
self.imap[iff],
#self.istep[iff],
self.time[iff],
) = self.read_fits_ifu(self.filenames[iff])
# Switch to a FORTRAN array
self.SI = np.asfortranarray(self.SI)
self.SQ = np.asfortranarray(self.SQ)
self.SU = np.asfortranarray(self.SU)
self.SV = np.asfortranarray(self.SV)
self.ICONT = np.asfortranarray(self.ICONT)
# Truncate unfullfilled maps
self.truncate_maps_ifu()
def truncate_maps(self):
"""Truncate unfullfilled maps"""
......@@ -1839,8 +1979,17 @@ class VFISVpackage:
else:
self.header['NSTEPS'] = max(self.istep[self.imap == self.header['NMAPS']]).astype("int")
def truncate_maps_ifu(self):
"""Truncate unfullfilled maps"""
def read_fits(self, filename):
# Change nmaps to only the number of maps that were fullfilled
self.header["NMAPS"] = max(self.imap).astype("int")
# introduce the RMAPS keyword for realized maps
self.header["RMAPS"] = self.header["NMAPS"]
def read_fits_slit(self, filename):
"""
Read the fits files into numpy array
......@@ -1878,6 +2027,45 @@ class VFISVpackage:
np.datetime64(hdul[0].header["DATE-OBS"]),
)
def read_fits_ifu(self, filename):
"""
Read the fits files into numpy array
Parameters
----------
filename : str
the fits filename to read
Returns
-------
Tuple of stokes vector and the imap and istep
"""
# Read the L1 split files and get the stokes vector and continuum intensity
# also get the imap and istep and datetime of the file
with fits.open(filename, memmap=True) as hdul:
data = hdul[0].data
si = data[0, self.pix_ini : self.pix_end + 1, :, :].copy().astype(np.float32)
sq = data[1, self.pix_ini : self.pix_end + 1, :, :].copy().astype(np.float32)
su = data[2, self.pix_ini : self.pix_end + 1, :, :].copy().astype(np.float32)
sv = data[3, self.pix_ini : self.pix_end + 1, :, :].copy().astype(np.float32)
ic = data[0, self.pix_cont, :, :].copy().astype(np.float32)
return (
si.T,
sq.T,
su.T,
sv.T,
ic.T,
hdul[0].header["IMAP"],
#hdul[0].header["ISTEP"],
np.datetime64(hdul[0].header["DATE-OBS"]),
)
def check_fits(self, header):
"""
Check's if the fits file is properly formatted.
......@@ -1922,7 +2110,7 @@ def main(path, id, out, numproc, line, width, preview, errors, diagnose, log):
if rank == 0:
# Starting Inversion
print_message(f"\nSDC:GRIS Inversion Pipeline v{get_version()}\n")
print_message(f"\nSDC:GRIS Inversion Pipeline: v{get_version()}\n")
# Check if the number of processor is passed via the mpirun -n flag or numproc
if not numproc:
......
Markdown is supported
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