Module themisasi.projections
Source code
import xarray
from pathlib import Path
import pymap3d as pm
import logging
import numpy as np
from matplotlib.pyplot import figure, draw, pause
from .plots import pcolormesh_nan, overlayrowcol
def asi_projection(dat: xarray.Dataset,
projalt_m: float = None,
min_el: float = 10.,
ofn: Path = None):
"""
plots ASI projected to altitude
* dat: Xarray containing image stack and metadata (az, el, lat, lon)
* projalt_m: projection altitude in meters
* min_el: minimum elevation angle (degrees). Data near the horizon is poorly calibrated (large angular error).
* ofn: filename to write of plot (optional)
"""
if projalt_m is None:
logging.error('projection altitude must be specified')
return
if dat['imgs'].shape[0] == 0:
return
if ofn:
ofn = Path(ofn).expanduser()
odir = ofn.parent
# %% censor pixels near the horizon with large calibration error do to poor skymap fits
badpix = dat['el'] < min_el
az = dat['az'].values
el = dat['el'].values
az[badpix] = np.nan
el[badpix] = np.nan
# %% coordinate transformation, let us know if error occurs
slant_range = projalt_m / np.sin(np.radians(el))
lat, lon, alt = pm.aer2geodetic(az, el, slant_range,
dat.lat.item(), dat.lon.item(), dat.alt_m.item())
# %% plots
fg = figure()
ax = fg.gca()
hi = pcolormesh_nan(lon, lat, dat['imgs'][0], cmap='gray', axis=ax) # priming
ttxt = f'Themis ASI {dat.site} projected to altitude {projalt_m/1e3} km\n' # FOV vs. HST0,HST1: green,red '
ht = ax.set_title(ttxt, color='g')
ax.set_xlabel('longitude')
ax.set_ylabel('latitude')
ax.autoscale(True, tight=True)
ax.grid(False)
# %% plot narrow FOV outline
if 'imgs2' in dat:
overlayrowcol(ax, dat.rows, dat.cols)
# %% play video
try:
for im in dat['imgs']:
ts = im.time.values.astype(str)[:-6]
hi.set_array(im.values.ravel()) # for pcolormesh
ht.set_text(ttxt + ts)
draw()
pause(0.01)
if ofn:
fn = odir / (ofn.stem + ts + ofn.suffix)
print('saving', fn, end='\r')
fg.savefig(fn, bbox_inches='tight', facecolor='k')
except KeyboardInterrupt:
return
def asi_radec(dat: xarray.Dataset,
min_el: float = 10.,
ofn: Path = None):
"""
plots ASI projected to altitude
* min_el: minimum elevation angle (degrees). Data near the horizon is poorly calibrated (large angular error).
* ofn: filename to write of plot (optional)
"""
if dat['imgs'].shape[0] == 0:
return
if ofn:
ofn = Path(ofn).expanduser()
az = dat.az.values
el = dat.el.values
az.setflags(write=True)
el.setflags(write=True)
bad = el < min_el # low to horizon, calibration is very bad
az[bad] = np.nan
el[bad] = np.nan
ra, dec = pm.azel2radec(az, el, dat.lat, dat.lon, dat.time)
fg = figure()
ax = fg.gca()
pcolormesh_nan(ra, dec, dat['imgs'].values[0], cmap='gray', axis=ax) # priming
ttxt = f'Themis ASI {dat.site} \n'
ax.set_title(ttxt, color='g')
ax.set_xlabel('right ascension [deg]')
ax.set_ylabel('declination [deg]')
ax.autoscale(True, tight=True)
ax.grid(False)
Functions
def asi_projection(dat, projalt_m=None, min_el=10.0, ofn=None)
-
plots ASI projected to altitude
- dat: Xarray containing image stack and metadata (az, el, lat, lon)
- projalt_m: projection altitude in meters
- min_el: minimum elevation angle (degrees). Data near the horizon is poorly calibrated (large angular error).
- ofn: filename to write of plot (optional)
Source code
def asi_projection(dat: xarray.Dataset, projalt_m: float = None, min_el: float = 10., ofn: Path = None): """ plots ASI projected to altitude * dat: Xarray containing image stack and metadata (az, el, lat, lon) * projalt_m: projection altitude in meters * min_el: minimum elevation angle (degrees). Data near the horizon is poorly calibrated (large angular error). * ofn: filename to write of plot (optional) """ if projalt_m is None: logging.error('projection altitude must be specified') return if dat['imgs'].shape[0] == 0: return if ofn: ofn = Path(ofn).expanduser() odir = ofn.parent # %% censor pixels near the horizon with large calibration error do to poor skymap fits badpix = dat['el'] < min_el az = dat['az'].values el = dat['el'].values az[badpix] = np.nan el[badpix] = np.nan # %% coordinate transformation, let us know if error occurs slant_range = projalt_m / np.sin(np.radians(el)) lat, lon, alt = pm.aer2geodetic(az, el, slant_range, dat.lat.item(), dat.lon.item(), dat.alt_m.item()) # %% plots fg = figure() ax = fg.gca() hi = pcolormesh_nan(lon, lat, dat['imgs'][0], cmap='gray', axis=ax) # priming ttxt = f'Themis ASI {dat.site} projected to altitude {projalt_m/1e3} km\n' # FOV vs. HST0,HST1: green,red ' ht = ax.set_title(ttxt, color='g') ax.set_xlabel('longitude') ax.set_ylabel('latitude') ax.autoscale(True, tight=True) ax.grid(False) # %% plot narrow FOV outline if 'imgs2' in dat: overlayrowcol(ax, dat.rows, dat.cols) # %% play video try: for im in dat['imgs']: ts = im.time.values.astype(str)[:-6] hi.set_array(im.values.ravel()) # for pcolormesh ht.set_text(ttxt + ts) draw() pause(0.01) if ofn: fn = odir / (ofn.stem + ts + ofn.suffix) print('saving', fn, end='\r') fg.savefig(fn, bbox_inches='tight', facecolor='k') except KeyboardInterrupt: return
def asi_radec(dat, min_el=10.0, ofn=None)
-
plots ASI projected to altitude
- min_el: minimum elevation angle (degrees). Data near the horizon is poorly calibrated (large angular error).
- ofn: filename to write of plot (optional)
Source code
def asi_radec(dat: xarray.Dataset, min_el: float = 10., ofn: Path = None): """ plots ASI projected to altitude * min_el: minimum elevation angle (degrees). Data near the horizon is poorly calibrated (large angular error). * ofn: filename to write of plot (optional) """ if dat['imgs'].shape[0] == 0: return if ofn: ofn = Path(ofn).expanduser() az = dat.az.values el = dat.el.values az.setflags(write=True) el.setflags(write=True) bad = el < min_el # low to horizon, calibration is very bad az[bad] = np.nan el[bad] = np.nan ra, dec = pm.azel2radec(az, el, dat.lat, dat.lon, dat.time) fg = figure() ax = fg.gca() pcolormesh_nan(ra, dec, dat['imgs'].values[0], cmap='gray', axis=ax) # priming ttxt = f'Themis ASI {dat.site} \n' ax.set_title(ttxt, color='g') ax.set_xlabel('right ascension [deg]') ax.set_ylabel('declination [deg]') ax.autoscale(True, tight=True) ax.grid(False)