Commit 47f2e23d authored by Vigeesh Gangadharan's avatar Vigeesh Gangadharan
Browse files

Merge branch 'develop' into 'master'

new previews

Closes #33

See merge request sdc/grisinv!16
parents b97ac9a7 93fce334
Pipeline #1429 passed with stage
in 2 minutes and 27 seconds
......@@ -15,15 +15,16 @@ import numpy as np
from datetime import date, datetime
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.gridspec as gridspec
import matplotlib.animation as animation
from scipy import interpolate
#supress fits obsgeo warnings
warnings.filterwarnings('ignore', category=FITSFixedWarning)
vfsiv_naming = {
vfisv_naming = {
"inte": {
"annotate": "S0+S1",
"unit": Unit.dimensionless_unscaled,
......@@ -54,7 +55,7 @@ vfsiv_naming = {
"incl": {
"annotate": "Inclination",
"unit": Unit.degree,
"filename": "inclincation",
"filename": "inclination",
"indices": [1],
"vmin": 0 * Unit.degree,
"vmax": 180 * Unit.degree,
......@@ -70,7 +71,38 @@ vfsiv_naming = {
"colormap": "BuPu",
},
}
vfsiv_errors_naming = {
vfisv_previews = {
"vlos": {
"annotate": "LOS Velocity",
"unit": Unit.cm / Unit.s,
"filename": "vlos",
"indices": [6],
"vmin": -3.0 * Unit.km / Unit.s,
"vmax": 3.0 * Unit.km / Unit.s,
"colormap": "RdBu",
},
"b_perp": {
"annotate": "|B$_{\perp}$|",
"unit": Unit.Gauss,
"filename": "bperp",
"indices": [None],
"vmin": 0. * Unit.Gauss,
"vmax": 3000. * Unit.Gauss,
"colormap": "summer",
},
"b_para": {
"annotate": "|B$_{\parallel}$|",
"unit": Unit.Gauss,
"filename": "bpara",
"indices": [None],
"vmin": 0. * Unit.Gauss,
"vmax": 3000. * Unit.Gauss,
"colormap": "summer",
}
}
vfisv_errors_naming = {
"chsq": {
"annotate": "Chi-sqr",
"unit": Unit.dimensionless_unscaled,
......@@ -150,26 +182,28 @@ def vfisv(
Returns
-------
data : dict of numpy arrays
inv_data : dict of numpy arrays
The dict contains the following variables
Continnum_I, vel_los, inclination, azimuth, field_strength
S0+S1,
vel_los, inclination, azimuth, field_strength and the corresponding uncertainities
chi_square.
stokes_data : dict of numpy arrays
The dict contains the following variables
si, sq,su,sv, si_fits, sq_fits, su_fits, sv_fits
header :
header : object
Astropy fits header
Examples
--------
from grisinv import vfisv
inv, header = vfisv('/dat/sdc/gris/20150920/level1_split/',5,15648.514,1.8,20)
inv, stokes, header = vfisv('/dat/sdc/gris/20150920/level1_split/',5,15648.514,1.8,20)
"""
vfisv_data = VFISVpackage(path, id, line, width)
# plt.plot(vfisv_data.WAVE, np.mean(vfisv_data.SI, axis=(0, 1)))
# plt.show()
# sys.exit()
# Spawn jobs to the grisinv:master+slaves
pycomm = MPI.COMM_SELF.Spawn(get_executable(), maxprocs=numproc)
......@@ -179,7 +213,7 @@ def vfisv(
for irank in range(numproc):
common_comm.Send([vfisv_data.line_props, MPI.FLOAT], dest=irank + 1, tag=5550)
common_comm.Send([vfisv_data.map_props, MPI.INT], dest=irank + 1, tag=5551)
common_comm.Send([vfisv_data.WAVE, MPI.FLOAT], dest=irank + 1, tag=7779)
common_comm.Send([vfisv_data.WAVE, MPI.FLOAT], dest=irank + 1, tag=5552)
# May cause issue: large buffers may create problem
for irank in range(numproc):
......@@ -189,9 +223,9 @@ def vfisv(
common_comm.Send([vfisv_data.SV, MPI.FLOAT], dest=irank + 1, tag=7773)
common_comm.Send([vfisv_data.ICONT, MPI.FLOAT], dest=irank + 1, tag=7774)
# receive the number of iterations from vfsiv:master
# receive the number of iterations from vfisv:master
niter = np.zeros(shape=(1), dtype="int32")
common_comm.Recv([niter, MPI.INT], source=1, tag=8888)
common_comm.Recv([niter, MPI.INT], source=1, tag=8880)
# progress through the iter
for i in tqdm(
......@@ -314,6 +348,7 @@ def get_istep_list(vfisv_params):
def get_result_shape(vfisv_params):
"""Get the shape of the data"""
nmaps = vfisv_params.header["NMAPS"]
nsteps = vfisv_params.header["NSTEPS"]
nlength = vfisv_params.header["NAXIS2"]
......@@ -323,6 +358,7 @@ def get_result_shape(vfisv_params):
def get_stokes_shape(vfisv_params):
"""Get the shape of the data"""
nmaps = vfisv_params.header["NMAPS"]
nsteps = vfisv_params.header["NSTEPS"]
nwave = vfisv_params.NUMW
......@@ -352,7 +388,7 @@ def get_quantity(data, vfisv_params, index: int, units=Unit.dimensionless_unscal
Returns
-------
quantity: float
the array in the proper shape
the array in the proper shape and units
"""
# Get the shape for the output
......@@ -389,7 +425,7 @@ def get_stokes_fits(data, vfisv_params):
data : float
array from vfisv containing all the data
vfisv_params :
the vfisvpackage class instance that hold all the params
the vfisvpackage class instance that holds all the params
index : int
the index of the desired variable based on VFISV manual
units :
......@@ -423,43 +459,49 @@ def get_stokes_fits(data, vfisv_params):
return quantity
def process_vfisv_output(data, sifit, sqfit, sufit, svfit, vfisv_params, errors=None):
def process_vfisv_output(data, sifit, sqfit, sufit, svfit, vfisv_params):
"""
Parameters
----------
data :
data : float
array from vfisv containing all the data
sifit: float
array containing the Stokes I fits
sqfit: float
array containing the Stokes Q fits
sufit: float
array containing the Stokes U fits
svfit: float
array containing the Stokes V fits
vfisv_params :
the vfisvpackage class instance that holds all the params
Returns
-------
inv_data : dict of numpy arrays
The dict contains the following variables
S0+S1, vel_los, inclination, azimuth, field_strength
stokes_data : dict of numpy arrays
The dict contains the following variables
si, sq,su,sv, si_fits, sq_fits, su_fits, sv_fits
"""
inv_data = {}
# get the data based on the index value and attach units.
for iquant in vfsiv_naming:
quant = 0
for i_index in vfsiv_naming[iquant]["indices"]:
quant += get_quantity(
data, vfisv_params, i_index, units=vfsiv_naming[iquant]["unit"]
)
inv_data[iquant] = quant
inv_data = add_data_from_vfisv_index(inv_data, data, vfisv_params, vfisv_naming)
# add errors
for iquant in vfsiv_errors_naming:
quant = 0
for i_index in vfsiv_errors_naming[iquant]["indices"]:
quant += get_quantity(
data, vfisv_params, i_index, units=vfsiv_errors_naming[iquant]["unit"]
)
inv_data[iquant] = quant
inv_data = add_data_from_vfisv_index(inv_data, data, vfisv_params, vfisv_errors_naming)
# TODO: Clean this also the flip
# Add the stokes vectors
stokes_data = {}
stokes_data["Stokes_I"] = get_stokes_fits(vfisv_params.SI, vfisv_params)
stokes_data["Stokes_Q"] = get_stokes_fits(vfisv_params.SQ, vfisv_params)
stokes_data["Stokes_U"] = get_stokes_fits(vfisv_params.SU, vfisv_params)
......@@ -470,19 +512,19 @@ def process_vfisv_output(data, sifit, sqfit, sufit, svfit, vfisv_params, errors=
stokes_data["Fitted_Stokes_U"] = get_stokes_fits(sufit, vfisv_params)
stokes_data["Fitted_Stokes_V"] = get_stokes_fits(svfit, vfisv_params)
#
# stokes_data["Stokes_I"] = vfisv_params.SI.T[:, ::-1, :]
# stokes_data["Stokes_Q"] = vfisv_params.SQ.T[:, ::-1, :]
# stokes_data["Stokes_U"] = vfisv_params.SU.T[:, ::-1, :]
# stokes_data["Stokes_V"] = vfisv_params.SV.T[:, ::-1, :]
# stokes_data["Fitted_Stokes_I"] = sifit.T[:, ::-1, :]
# stokes_data["Fitted_Stokes_Q"] = sqfit.T[:, ::-1, :]
# stokes_data["Fitted_Stokes_U"] = sufit.T[:, ::-1, :]
# stokes_data["Fitted_Stokes_V"] = svfit.T[:, ::-1, :]
return inv_data, stokes_data
def add_data_from_vfisv_index(inv_data, data,vfisv_params,vfisv_dict):
"Add the data into array based on the VFISV indices"
for iquant in vfisv_dict:
quant = 0
for i_index in vfisv_dict[iquant]["indices"]:
quant += get_quantity(
data, vfisv_params, i_index, units=vfisv_dict[iquant]["unit"]
)
inv_data[iquant] = quant
return inv_data
def write_stokes_fits(stokes_data, header, diagnose):
......@@ -559,7 +601,7 @@ def write_fits(data, header, out, errors=None):
header["FILENAME"] = out.split("/")[-1]
hdu_list = []
for iquant in vfsiv_naming:
for iquant in vfisv_naming:
# intensity is the primary hdu, all the other quantities are image hdu
if iquant == "inte":
hdu = fits.PrimaryHDU(data=data[iquant].value, header=header)
......@@ -567,13 +609,13 @@ def write_fits(data, header, out, errors=None):
hdu = fits.ImageHDU(data=data[iquant].value, header=header)
# adding the unique extname
hdu.header["EXTNAME"] = extname + "_" + vfsiv_naming[iquant]["filename"]
hdu.header["EXTNAME"] = extname + "_" + vfisv_naming[iquant]["filename"]
# adding the description
hdu.header["BTYPE"] = vfsiv_naming[iquant]["annotate"]
hdu.header["BTYPE"] = vfisv_naming[iquant]["annotate"]
# adding the unit of the qauntity
hdu.header["BUNIT"] = str(vfsiv_naming[iquant]["unit"])
hdu.header["BUNIT"] = str(vfisv_naming[iquant]["unit"])
# append the hdu in the list
hdu_list.append(hdu)
......@@ -592,7 +634,7 @@ def write_fits(data, header, out, errors=None):
if errors != None:
hdu_list = []
for iquant in vfsiv_errors_naming:
for iquant in vfisv_errors_naming:
# intensity is the primary hdu, all the other quantities are image hdu
if iquant == "chsq":
hdu = fits.PrimaryHDU(data=data[iquant].value, header=header)
......@@ -601,14 +643,14 @@ def write_fits(data, header, out, errors=None):
# adding the unique extname
hdu.header["EXTNAME"] = (
extname + "_" + vfsiv_errors_naming[iquant]["filename"]
extname + "_" + vfisv_errors_naming[iquant]["filename"]
)
# adding the description
hdu.header["BTYPE"] = vfsiv_errors_naming[iquant]["annotate"]
hdu.header["BTYPE"] = vfisv_errors_naming[iquant]["annotate"]
# adding the unit of the qauntity
hdu.header["BUNIT"] = str(vfsiv_errors_naming[iquant]["unit"])
hdu.header["BUNIT"] = str(vfisv_errors_naming[iquant]["unit"])
# append the hdu in the list
hdu_list.append(hdu)
......@@ -649,7 +691,7 @@ def add_grisinv_info(hdr, vfisv_params):
hdr["WAVERPIX"] = 0
hdr.comments["WAVERPIX"] = "Wavelength reference pixel"
hdr["WAVERVAL"] = vfisv_params.WAVE[0] * 1000.0 # mAngstrom->Angstrom
hdr["WAVERVAL"] = vfisv_params.wave_start
hdr.comments["WAVERVAL"] = "Wavelength reference value"
hdr["WAVEDELT"] = vfisv_params.header["CDELT3"]
......@@ -833,14 +875,28 @@ def make_header(vfisv_params):
return hdr
def get_cbar_ax(ax):
divider = make_axes_locatable(ax)
return divider.append_axes("top", size="5%", pad=0.05)
def plot_image(
data, header, preview, cbar_title=None, dpi=100, vmin=None, vmax=None, cmap=None
):
"""
Plot preview images of the inversion data based on the header params
Parameters
----------
data : dict of numpy arrays
The dict contains the following variables
Continnum_I, vel_los, inclination, azimuth, field_strength
header: object
Astropy fits header
preview: str
filename prefix for the preview images
Returns
-------
"""
# TODO: Clean this. make is less cluttered.
wcs = WCS(header)
......@@ -854,7 +910,10 @@ def plot_image(
aspect_ratio = data.shape[1] / data.shape[2]
shape_inches = (max([8 / aspect_ratio + 0.5, 3]) + 1.0, 8)
fig = plt.figure(figsize=shape_inches, dpi=dpi)
ax = plt.subplot(projection=wcs, slices=("x", "y", 0))
spec = gridspec.GridSpec(ncols=3, nrows=2, figure=fig, width_ratios=[data.shape[2],0,5.])
ax = fig.add_subplot(spec[0:2, 0], projection=wcs, slices=["x", "y", 0])
cax = fig.add_subplot(spec[1, 2])
data_units = data.unit
if header["NMAPS"] == 1:
if vmax == None:
unit_data = data[0, :, :].value
......@@ -862,6 +921,7 @@ def plot_image(
unit_vmax = vmax
else:
unit_data = data[0, :, :].to(vmax.unit).value
data_units = data[0, :, :].to(vmax.unit).unit
unit_vmin = vmin.value
unit_vmax = vmax.value
......@@ -875,14 +935,21 @@ def plot_image(
cb = ax.imshow(
unit_data, origin="lower", vmin=unit_vmin, vmax=unit_vmax, cmap=cmap
)
ax.grid(color="dimgray", ls="dashed")
ax.set_xlabel("Heliprojective Longitude (arcsec)")
ax.set_ylabel("Heliprojective Latitude (arcsec)")
ax.set_title(title)
cbar = plt.colorbar(cb)
cbar.ax.set_ylabel(cbar_title)
plt.tight_layout()
fig.subplots_adjust(bottom=0.1, top=0.9)
#ax.set_title(title)
fig.suptitle(title)
# set lon_lat axis
set_axis_lon_lat(ax)
# additonal axis
set_additional_axis(ax, header)
cbar = plt.colorbar(cb, cax=cax)
cbar.ax.set_ylabel(cbar_title+ f" ({data_units})")
fig.canvas.draw()
fig.tight_layout()
fig.subplots_adjust(bottom=0.1, top=0.85, left=1/shape_inches[0], right=1-1/shape_inches[0])
plt.savefig(preview) # .replace('.',f'_{imap:03d}.'))
# plt.clf()
else:
......@@ -895,6 +962,7 @@ def plot_image(
unit_vmax = vmax
else:
unit_data = data[imap, :, :].to(vmax.unit).value
data_units = data[0, :, :].to(vmax.unit).unit
unit_vmin = vmin.value
unit_vmax = vmax.value
......@@ -918,37 +986,102 @@ def plot_image(
ax.imshow(
unit_data, origin="lower", vmin=unit_vmin, vmax=unit_vmax, cmap=cmap
)
ax.grid(color="dimgray", ls="dashed")
ax.set_xlabel("Heliprojective Longitude (arcsec)")
ax.set_ylabel("Heliprojective Latitude (arcsec)")
cbar = plt.colorbar(cb)
cbar.ax.set_ylabel(cbar_title)
plt.tight_layout()
fig.subplots_adjust(bottom=0.1, top=0.9)
ax.set_title(title)
# set lon_lat axis
set_axis_lon_lat(ax)
# additonal axis
set_additional_axis(ax,header)
cbar = plt.colorbar(cb, cax=cax)
cbar.ax.set_ylabel(cbar_title+ f" ({data_units})")
fig.canvas.draw()
fig.tight_layout()
fig.subplots_adjust(bottom=0.1, top=0.85, left=1/shape_inches[0],right=1-1/shape_inches[0])
fig.suptitle(title)
#ax.set_title(title)
ims.append([cb])
ani = animation.ArtistAnimation(fig, ims, interval=1000.0, blit=True)
ani.save(preview.replace(".png", ".gif"))
if ".png" in preview:
ani.save(preview.replace(".png", ".gif"))
elif ".jpeg" in preview:
ani.save(preview.replace(".jpeg", ".gif"))
elif ".jpg" in preview:
ani.save(preview.replace(".jpeg", ".gif"))
else:
ani.save(preview+".gif")
return None
def set_axis_lon_lat(ax):
"""Format the Longitude and Latitide axes"""
lon = ax.coords['hpln']
lat = ax.coords['hplt']
lon.set_axislabel("HPLN-TAN")
lat.set_axislabel("HPLT-TAN")
lon.set_ticks(spacing=5 * Unit.arcsec, color='dimgray')
lat.set_ticks(spacing=5 * Unit.arcsec, color='dimgray')
lon.grid(color='dimgray', alpha=0.5, linestyle='dashed')
lat.grid(color='dimgray', alpha=0.5, linestyle='dashed')
lon.tick_params(
axis='both',
which='both',
bottom=True,
top=False,
right=False,
left=False,
labelbottom=True)
lat.tick_params(
axis='both',
which='both',
bottom=False,
top=False,
right=False,
left=True,
labelleft=True)
return None
def set_additional_axis(ax,header):
"""Format the x/y axes"""
xdeltarcsec = ((abs(header['CDELT1'])) * Unit.degree).to(Unit.arcsecond).value
ydeltarcsec = ((abs(header['CDELT2'])) * Unit.degree).to(Unit.arcsecond).value
secax_x = ax.secondary_xaxis('top', functions=(lambda x: x * xdeltarcsec, lambda x: x / xdeltarcsec))
secax_y = ax.secondary_yaxis('right', functions=(lambda x: x * ydeltarcsec, lambda x: x / ydeltarcsec))
secax_x.set_xlabel("x (arcsec)")
secax_y.set_ylabel("y (arcsec)")
def create_preview(data, header, preview):
"""Create preview image"""
for iquant in vfisv_previews:
if iquant=='b_para':
quant_data = np.abs(data['bmag']*np.cos(data['incl']))
elif iquant=='b_perp':
quant_data = np.abs(data['bmag']*np.sin(data['incl']))
else:
quant_data = data[iquant]
for iquant in vfsiv_naming:
quant_data = data[iquant]
# velocity calibration for preview images only
if iquant == "vlos":
quant_data -= np.mean(quant_data)
#setting a magnetic field mask of 200 Guass.
mag_mask = data["bmag"].value < 200.
quant_data -= np.mean(quant_data[mag_mask])
plot_image(
quant_data,
header,
preview.replace(".", f"_{vfsiv_naming[iquant]['filename']}."),
cbar_title=vfsiv_naming[iquant]["annotate"]
+ f" ({vfsiv_naming[iquant]['unit']})",
vmin=vfsiv_naming[iquant]["vmin"],
vmax=vfsiv_naming[iquant]["vmax"],
cmap=vfsiv_naming[iquant]["colormap"],
preview.replace(".", f"_{vfisv_previews[iquant]['filename']}_preview."),
cbar_title=vfisv_previews[iquant]["annotate"],
vmin=vfisv_previews[iquant]["vmin"],
vmax=vfisv_previews[iquant]["vmax"],
cmap=vfisv_previews[iquant]["colormap"],
)
return None
......@@ -1054,6 +1187,7 @@ class VFISVpackage:
self.date_beg = None
self.date_end = None
self.fits_wcs = None
self.wave_start = None
# Populate all parameters
self.get_filenames()
......@@ -1181,6 +1315,8 @@ class VFISVpackage:
.value
)
self.wave_start= wave_array[0]
if (wave_array[0] < self.line) and (wave_array[-1] > self.line):
wave_func = interpolate.interp1d(wave_array, np.arange(len(wave_array)))
self.pix_lab_wav = wave_func(self.line) # findindex(wave_array, self.line)
......@@ -1250,6 +1386,10 @@ class VFISVpackage:
self.ICONT = np.asfortranarray(self.ICONT)
# Change nmaps to reflect the fullfilled maps
self.header['NMAPS'] = max(self.imap).astype("int")
def read_fits(self, filename):
"""
Read the fits files into numpy array
......
......@@ -7,7 +7,7 @@ F77F1 = -c -Ofast -std=legacy #-fcheck=all -g -Wall
F77F2 = -c -Ofast -std=legacy #-fcheck=all -g -Wall
#Pipeline under conda
LIBS = -L $CONDA_PREFIX/lib -llapack
LIBS = -L$(CONDA_PREFIX)/lib -llapack
OBJ = hold.o \
time_param.o \
......
......@@ -115,7 +115,7 @@ PROGRAM VFISV_SPEC
! Read spectral line data, Keep for allocation.
CALL READ_SPEC_LINE(NTASKS)
! Initialize Wavelength grid
CALL MPI_RECV(WAVE,NUMW,MPI_REAL,0,7779,intracomm,STATUS,IERR)
CALL MPI_RECV(WAVE,NUMW,MPI_REAL,0,5552,intracomm,STATUS,IERR)
! Keep for allocation
CALL READ_SPEC_MAP(NTASKS)
......@@ -156,7 +156,7 @@ PROGRAM VFISV_SPEC
ITERR=INT(DBLE(LX)/DBLE(MYTASK-2))
!
!Send the iter to the grand master.
IF (NTASKS == 1) CALL MPI_SEND(ITERR,1,MPI_INTEGER,0,8888,intracomm,IERR)
IF (NTASKS == 1) CALL MPI_SEND(ITERR,1,MPI_INTEGER,0,8880,intracomm,IERR)
DO J=1,ITERR+1
IF (J <= ITERR) THEN
......
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