Commit bc2766a1 authored by Vigeesh Gangadharan's avatar Vigeesh Gangadharan
Browse files

Merge branch 'develop' into 'master'

Changes to preview images, temporary crpix fix.

See merge request sdc/grisinv!20
parents 6b5f5432 3cc78864
Pipeline #1477 passed with stage
in 2 minutes and 5 seconds
......@@ -20,7 +20,6 @@ import matplotlib.animation as animation
from scipy import interpolate
#supress fits obsgeo warnings
warnings.filterwarnings('ignore', category=FITSFixedWarning)
......@@ -80,25 +79,34 @@ vfisv_previews = {
"indices": [6],
"vmin": -3.0 * Unit.km / Unit.s,
"vmax": 3.0 * Unit.km / Unit.s,
"colormap": "RdBu",
"colormap": "RdBu_r",
"gridcolor" : "black",
"label_fmt" : "%1.1f",
"cticks" : ([-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0])
},
"b_perp": {
"annotate": "|B$_{\perp}$|",
"unit": Unit.Gauss,
"filename": "bperp",
"indices": [None],
"vmin": 0. * Unit.Gauss,
"vmax": 3000. * Unit.Gauss,
"colormap": "summer",
"vmin": 0. * Unit.kiloGauss,
"vmax": 3. * Unit.kiloGauss,
"colormap": "viridis",
"gridcolor" : "white",
"label_fmt" : "%1.1f",
"cticks" : ([0.0,0.5,1.0,1.5,2.0,2.5,3.0])
},
"b_para": {
"annotate": "|B$_{\parallel}$|",
"unit": Unit.Gauss,
"filename": "bpara",
"indices": [None],
"vmin": 0. * Unit.Gauss,
"vmax": 3000. * Unit.Gauss,
"colormap": "summer",
"vmin": 0. * Unit.kiloGauss,
"vmax": 3. * Unit.kiloGauss,
"colormap": "viridis",
"gridcolor" : "white",
"label_fmt" : "%1.1f",
"cticks" : ([0.0,0.5,1.0,1.5,2.0,2.5,3.0])
}
}
......@@ -688,7 +696,7 @@ def add_grisinv_info(hdr, vfisv_params):
hdr["WAVE_STR"] = vfisv_params.line_name
hdr.comments["WAVE_STR"] = "Spectral line name"
hdr["WAVERPIX"] = 0
hdr["WAVERPIX"] = 1
hdr.comments["WAVERPIX"] = "Wavelength reference pixel"
hdr["WAVERVAL"] = vfisv_params.wave_start
......@@ -715,14 +723,13 @@ def add_grisinv_info(hdr, vfisv_params):
hdr["PRPROC1"] = "VFISV"
hdr.comments["PRPROC1"] = "Name of procedure used"
hdr[
"PRPARA1"
] = f'PATH="{vfisv_params.path}",ID={vfisv_params.id},LINE={vfisv_params.line},WIDTH={vfisv_params.width}'
prpara1 = f'PATH="{vfisv_params.path}",ID={vfisv_params.id},LINE={vfisv_params.line},WIDTH={vfisv_params.width}'
hdr["PRPARA1"] = prpara1
hdr.comments["PRPARA1"] = "List of parameters/options for PRPROCn"
# hdr['PRBRA1'] =
# hdr.comments['PRBRA1'] =
#
# hdr['PRVER1A'] =
# hdr.comments['PRVER1A'] =
......@@ -876,7 +883,7 @@ def make_header(vfisv_params):
def plot_image(
data, header, preview, cbar_title=None, dpi=100, vmin=None, vmax=None, cmap=None
data, header, preview, cbar_title=None, dpi=100, vmin=None, vmax=None, cmap=None, grid=None, cfmt=None,cticks=None
):
"""
Plot preview images of the inversion data based on the header params
......@@ -898,7 +905,7 @@ def plot_image(
"""
# TODO: Clean this. make is less cluttered.
# TODO: Clean this. make it less cluttered.
wcs = WCS(header)
title = (
......@@ -909,28 +916,22 @@ 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)
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
fig = plt.figure(figsize=shape_inches, dpi=dpi, constrained_layout=True)
spec = gridspec.GridSpec(nrows=2, ncols=2, figure=fig, width_ratios=[data.shape[2],5.],wspace=0.1,hspace=0.)
ax = fig.add_subplot(spec[:, 0], projection=wcs, slices=["x", "y", 0])
cax = fig.add_subplot(spec[1, 1])
if header["NMAPS"] == 1:
if vmax == None:
unit_data = data[0, :, :].value
unit_vmin = vmin
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
unit_data, data_units, unit_vmin, unit_vmax = convert_units(data,vmin,vmax,0)
# a quick correction for plots of field
unit_vmax = (
np.max(unit_data)
if (unit_vmax != None) and (unit_vmax > np.max(unit_data))
else unit_vmax
)
if data[0, ...].unit == Unit.Gauss:
unit_vmax = (
np.max(unit_data)
if (unit_vmax != None) and (unit_vmax > np.max(unit_data))
else unit_vmax
)
cb = ax.imshow(
unit_data, origin="lower", vmin=unit_vmin, vmax=unit_vmax, cmap=cmap
......@@ -940,34 +941,23 @@ def plot_image(
fig.suptitle(title)
# set lon_lat axis
set_axis_lon_lat(ax)
set_axis_lon_lat(ax,grid=grid)
# additonal axis
set_additional_axis(ax, header)
cbar = plt.colorbar(cb, cax=cax)
cbar = plt.colorbar(cb, cax=cax,format=cfmt)
cbar.set_ticks(cticks)
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()
plt.savefig(preview,bbox_inches='tight')
else:
ims = []
for imap in range(header["NMAPS"]):
# TODO: Clean this
if vmax == None:
unit_data = data[imap, :, :].value
unit_vmin = vmin
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
unit_data, data_units, unit_vmin, unit_vmax = convert_units(data,vmin,vmax,imap)
# a quick correction for plots of field
if data[imap, :, :].unit == Unit.Gauss:
if data[imap, ...].unit == Unit.Gauss:
unit_vmax = (
np.max(unit_data)
if (unit_vmax != None) and (unit_vmax > np.max(unit_data))
......@@ -983,29 +973,23 @@ def plot_image(
animated=True,
)
if imap == 0:
ax.imshow(
unit_data, origin="lower", vmin=unit_vmin, vmax=unit_vmax, cmap=cmap
)
# set lon_lat axis
set_axis_lon_lat(ax)
set_axis_lon_lat(ax,grid=grid)
# additonal axis
set_additional_axis(ax,header)
cbar = plt.colorbar(cb, cax=cax)
cbar = plt.colorbar(cb, cax=cax,format=cfmt)
cbar.set_ticks(cticks)
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 = animation.ArtistAnimation(fig, ims, interval=500.0, blit=True)
if ".png" in preview:
ani.save(preview.replace(".png", ".gif"))
elif ".jpeg" in preview:
......@@ -1017,7 +1001,26 @@ def plot_image(
return None
def set_axis_lon_lat(ax):
def convert_units(data,vmin,vmax,imap):
'''Convert the units based on vmin/vmax'''
if vmax == None:
# if vmax is not provided
unit_data = data[imap, ...].value
unit_vmin = vmin
unit_vmax = vmax
data_units = data.unit
else:
# if vmax is not provided with a unit.
unit_data = data[imap, ...].to(vmax.unit).value
data_units = data[imap, ...].to(vmax.unit).unit
unit_vmin = vmin.value
unit_vmax = vmax.value
return unit_data, data_units, unit_vmin, unit_vmax
def set_axis_lon_lat(ax,grid=None):
"""Format the Longitude and Latitide axes"""
lon = ax.coords['hpln']
lat = ax.coords['hplt']
......@@ -1025,8 +1028,9 @@ def set_axis_lon_lat(ax):
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')
if grid:
lon.grid(color=grid, alpha=0.5, linestyle='dashed')
lat.grid(color=grid, alpha=0.5, linestyle='dashed')
lon.tick_params(
axis='both',
......@@ -1046,17 +1050,31 @@ def set_axis_lon_lat(ax):
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
xdeltarcsec = attach_unit(abs(header['CDELT1']),Unit.degree,convert=Unit.arcsecond,value=True)
ydeltarcsec = attach_unit(abs(header['CDELT2']), Unit.degree, convert=Unit.arcsecond, value=True)
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 attach_unit(data,units,convert=None,value=False):
'''Attach the units to the variable and if given convert'''
# attach units
data_units=data*units
#convert to the requested units
if convert:
data_units = data_units.to(convert)
# if only the value is required
if value:
data_units = data_units.value
return data_units
def create_preview(data, header, preview):
"""Create preview image"""
......@@ -1070,8 +1088,8 @@ def create_preview(data, header, preview):
# velocity calibration for preview images only
if iquant == "vlos":
#setting a magnetic field mask of 200 Guass.
mag_mask = data["bmag"].value < 200.
#setting a magnetic field mask of 1000 Guass.
mag_mask = data["bmag"].value < 1000.
quant_data -= np.mean(quant_data[mag_mask])
plot_image(
......@@ -1082,6 +1100,9 @@ def create_preview(data, header, preview):
vmin=vfisv_previews[iquant]["vmin"],
vmax=vfisv_previews[iquant]["vmax"],
cmap=vfisv_previews[iquant]["colormap"],
grid=vfisv_previews[iquant]["gridcolor"],
cfmt=vfisv_previews[iquant]["label_fmt"],
cticks=vfisv_previews[iquant]["cticks"]
)
return None
......@@ -1195,7 +1216,7 @@ class VFISVpackage:
self.get_data()
self.get_map_props()
self.get_line_props()
self.get_wave()
self.get_rel_wave()
self.nmaps = self.header["NMAPS"]
self.nsteps = self.header["NSTEPS"] # does the number of steps per map change
......@@ -1203,17 +1224,25 @@ class VFISVpackage:
self.LY = self.header["NAXIS2"]
self.LZ = 16 # VFISV data structure
def get_line_list(self):
"""Get the list of available lines"""
line_list = (
pkgutil.get_data("grisinv", os.path.join("data_files", "line_names.csv"))
.decode()
.split("\n")
)
return line_list
def check_spectral_line(self):
"""
Check if the spectral line is in the list of lines
If the spectral line exists, then get the geff value.
"""
line_list = (
pkgutil.get_data("grisinv", os.path.join("data_files", "line_names.csv"))
.decode()
.split("\n")
)
line_list = self.get_line_list()
for iline in line_list:
# Todo: change this to numpy float comapre.
if self.line == float(iline.split(",")[1]):
self.line_name = iline.split(",")[0]
self.geff = iline.split(",")[3]
......@@ -1258,14 +1287,41 @@ class VFISVpackage:
dtype="int32",
)
def get_wave(self):
"""Get the wavlength array from the header"""
stepww = self.header["CDELT3"] * 1000.0 # in mAngstrom
self.NUMW = np.array(self.pix_end - self.pix_ini + 1, dtype="int32")
self.WAVE = (
np.arange(self.NUMW, dtype="float32") * stepww
- (self.pix_lab_wav - self.pix_ini) * stepww
)
def get_wavelength_array(self,convert=None):
"""Get the wavelength array from the WCS keywords"""
# Get the WCS of the wavelength axis
wcs_wavelength = WCS(self.header, naxis=[3])
# Create pixel array of size wavelength array
wavelength_pix = np.arange(self.header["NAXIS3"])
# Convert from pixels to world
# TODO: Change the origin to 0
wave_array = wcs_wavelength.wcs_pix2world(wavelength_pix, 1)[0]
# Convert the wavelength array into the required uints and return value
if convert:
wave_array *= Unit.m
return wave_array.to(convert).value
else:
return wave_array
def get_rel_wave(self):
"""Get the wavelength relative to the line for VFISV"""
# Get the wavelength array from WCS
wave_array = self.get_wavelength_array(convert=Unit.Angstrom)
# Compute the array relative to the spectral line
wave_array -= self.line
# VFISV needs this array in mAngstrom
wave_array *= 1000.0 # in mAngstrom
# Set the wavelength and num_wavelength for VFISV
self.WAVE = wave_array[self.pix_ini:self.pix_end+1].astype("float32")
self.NUMW = np.array(len(self.WAVE), dtype="int32")
def get_filenames(self):
"""Get filenames from the data and observation id"""
......@@ -1299,32 +1355,31 @@ class VFISVpackage:
self.date_end = datetime.fromisoformat(header_end["DATE-BEG"])
def get_pixel_bounds(self):
"""Get the +/- width A bounds centered at the line"""
# wave_array = (
# self.header["CRVAL3"]
# + np.arange(self.header["NAXIS3"]) * self.header["CDELT3"]
# )
wave_array = (
(
WCS(self.header, naxis=[3]).wcs_pix2world(
np.arange(self.header["NAXIS3"]), 0
)[0]
* Unit.m
)
.to(Unit.Angstrom)
.value
)
"""Get the +/- width Angstrom bounds centered at the line"""
wave_array = self.get_wavelength_array(convert=Unit.Angstrom)
# Check if the line is in limit
if wave_array[0] <= self.line <= wave_array[-1]:
# Get wavelength relative to the spectral line
wave_array_rel = wave_array - self.line
# Interpolate to find the float-index of the spectral line
wave_func = interpolate.interp1d(wave_array_rel, np.arange(len(wave_array_rel)))
self.pix_lab_wav = wave_func(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)
#get the ini and end pix
self.pix_ini, self.pix_end = (
int(self.pix_lab_wav - self.width // self.header["CDELT3"]),
int(self.pix_lab_wav + self.width // self.header["CDELT3"]),
np.round(self.pix_lab_wav - self.width / self.header["CDELT3"]).astype('int'),
np.round(self.pix_lab_wav + self.width / self.header["CDELT3"]).astype('int'),
)
# Set the pix continnum to pix_end
self.pix_cont = self.pix_end
# Set the start of the wavelenght array for
self.wave_start = wave_array[self.pix_ini]
else:
......@@ -1371,10 +1426,10 @@ class VFISVpackage:
+ "{percentage:3.0f}% ({elapsed_s:3.1f}s)",
):
(
self.SI[iff, :, :],
self.SQ[iff, :, :],
self.SU[iff, :, :],
self.SV[iff, :, :],
self.SI[iff, ...],
self.SQ[iff, ...],
self.SU[iff, ...],
self.SV[iff, ...],
self.ICONT[iff, :],
self.imap[iff],
self.istep[iff],
......@@ -1409,10 +1464,10 @@ class VFISVpackage:
data = hdul[0].data
si = data[0, self.pix_ini - 1 : self.pix_end, :, 0].copy()
sq = data[1, self.pix_ini - 1 : self.pix_end, :, 0].copy()
su = data[2, self.pix_ini - 1 : self.pix_end, :, 0].copy()
sv = data[3, self.pix_ini - 1 : self.pix_end, :, 0].copy()
si = data[0, self.pix_ini : self.pix_end + 1, :, 0].copy()
sq = data[1, self.pix_ini : self.pix_end + 1, :, 0].copy()
su = data[2, self.pix_ini : self.pix_end + 1, :, 0].copy()
sv = data[3, self.pix_ini : self.pix_end + 1, :, 0].copy()
ic = data[0, self.pix_cont, :, 0].copy()
return (
......
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