From 093d43934d03ccd0e494a29328d7615ff6c19f67 Mon Sep 17 00:00:00 2001 From: Edward Molter Date: Mon, 2 Oct 2023 08:37:47 -0700 Subject: [PATCH 1/2] adding io support and some nav functionality --- pylanetary/navigation/core.py | 22 +++++++++---- pylanetary/utils/core.py | 61 ++++++++++++++++++++++++----------- pylanetary/utils/data/Io.yaml | 4 +-- 3 files changed, 61 insertions(+), 26 deletions(-) diff --git a/pylanetary/navigation/core.py b/pylanetary/navigation/core.py index bd6af07..12c5f46 100755 --- a/pylanetary/navigation/core.py +++ b/pylanetary/navigation/core.py @@ -256,7 +256,7 @@ def limb_darkening(mu, a, law='exp', mu0=None): https://www.astro.keele.ac.uk/jkt/codes/jktld.html ''' # Input check - if np.any(mu <0) or np.any(mu>1): + if np.any(mu <0.0) or np.any(mu>1.0): raise ValueError('Cosine of emission angle range [0,1]') mu = np.array(mu) #necessary so when floats are passed in, the line mu[bad] = 0 doesnt complain about indexing a float @@ -388,6 +388,10 @@ def __init__(self, ob_lon, ob_lat, pixscale_km, np_ang, req, rpol, center=(0.0, self.sun_n = surface_normal(self.lat_g, self.lon_w, self.sun_lon) self.mu0 = emission_angle(self.ob_lat, self.sun_n) + # solve small numerical issue where mu = 1.0 + epsilon + self.mu[(self.mu > 1.0)*(self.mu < 1.00001)] = 1.0 + self.mu0[(self.mu0 > 1.0)*(self.mu0 < 1.00001)] = 1.0 + avg_circumference = 2*np.pi*((self.req + self.rpol)/2.0) self.deg_per_px = self.pixscale_km * (1/avg_circumference) * 360 @@ -403,7 +407,7 @@ def __str__(self): return f'ModelEllipsoid instance; req={self.req}, rpol={self.rpol}' - def ldmodel(self, tb, a, beam=None, law='exp', mu0=None): + def ldmodel(self, tb, a, beam=None, law='exp', mu0=None, psf_mode='gaussian'): ''' Make a limb-darkened model disk convolved with the beam See docstring of limb_darkening() for options @@ -424,6 +428,8 @@ def ldmodel(self, tb, a, beam=None, law='exp', mu0=None): has no effect unless law=="minnaert". if self.mu0 undefined, law=="minnaert", and mu0=None, then code will fail + psf_mode : str, optional, default "gaussian" + mode to use for convolve_with_beam(), options "airy", "gaussian" ''' if mu0 is None: mu0=getattr(self, "mu0", None) @@ -436,7 +442,7 @@ def ldmodel(self, tb, a, beam=None, law='exp', mu0=None): if beam is None: return ldmodel else: - return convolve_with_beam(ldmodel, beam) + return convolve_with_beam(ldmodel, beam, mode=psf_mode) def zonalmodel(lats, lons, tbs, a=0.0): @@ -573,6 +579,8 @@ def colocate(self, mode = 'convolution', diagnostic_plot=True, save_plot=None, * type of limb darkening model to use :beam: float, tuple, or np.array, optional. default None. units pixels. see utils.convolve_with_beam + :psf_mode: str, optional. default "gaussian" + what beam shape to use, options "airy", "gaussian" :low_thresh: float, required. see documentation of skimage.feature.canny :high_thresh: float, required. @@ -594,6 +602,8 @@ def colocate(self, mode = 'convolution', diagnostic_plot=True, save_plot=None, * type of limb darkening model to use :beam: float, tuple, or np.array, optional. default None. units pixels. see utils.convolve_with_beam + :psf_mode: str, optional. default "gaussian" + what beam shape to use, options "airy", "gaussian" :err: float per-pixel error in input image @@ -628,17 +638,17 @@ def colocate(self, mode = 'convolution', diagnostic_plot=True, save_plot=None, * ----- * sometimes dxerr, dyerr give unrealistic or undefined behavior ''' - defaultKwargs={'err':None,'beam':None,'law':'exp'} + defaultKwargs={'err':None,'beam':None,'law':'exp', 'psf_mode':'gaussian'} kwargs = { **defaultKwargs, **kwargs } if (mode == 'convolution') or (mode == 'disk'): model = self.ldmodel(kwargs['tb'], kwargs['a'], law=kwargs['law']) - model = convolve_with_beam(model, kwargs['beam']) + model = convolve_with_beam(model, kwargs['beam'], mode=kwargs['psf_mode']) data_to_compare = self.data elif mode == 'canny': model_planet = self.ldmodel(kwargs['tb'], kwargs['a'], law=kwargs['law']) if kwargs['beam'] is not None: - model_planet = convolve_with_beam(model_planet, kwargs['beam']) + model_planet = convolve_with_beam(model_planet, kwargs['beam'], mode=kwargs['psf_mode']) edges = feature.canny(self.data, sigma=kwargs['sigma'], low_threshold = kwargs['low_thresh'], high_threshold = kwargs['high_thresh']) model = feature.canny(model_planet, sigma=kwargs['sigma'], low_threshold = kwargs['low_thresh'], high_threshold = kwargs['high_thresh']) diff --git a/pylanetary/utils/core.py b/pylanetary/utils/core.py index 3e4944f..7d6d592 100644 --- a/pylanetary/utils/core.py +++ b/pylanetary/utils/core.py @@ -1,16 +1,18 @@ # Licensed under a ??? style license - see LICENSE.rst import numpy as np +from datetime import timedelta +import importlib, yaml, importlib.resources +from scipy.interpolate import interp1d +from scipy import fftpack + from astropy import convolution from astropy.coordinates import Angle import astropy.units as u from astropy.units import Quantity -import importlib, yaml, importlib.resources -from scipy.interpolate import interp1d -from scipy import fftpack -import pylanetary.utils.data as body_info from astropy.time import Time from astroquery.jplhorizons import Horizons -from datetime import timedelta + +import pylanetary.utils.data as body_info ''' To do: @@ -363,7 +365,7 @@ def rebin(arr, z): return (1/z**2)*arr.reshape(shape).mean(-1).mean(1) -def convolve_with_beam(data, beam): +def convolve_with_beam(data, beam, mode='gaussian'): ''' Convolves input 2-D image with a Gaussian beam or an input PSF image @@ -371,34 +373,57 @@ def convolve_with_beam(data, beam): ---------- data : np.array, required. beam : float/int, 3-element array-like, np.array, or None, required. - if float/int, circular Gaussian beam assumed, and this sets the fwhm - [pixels]. - if 3-element array-like, those are (fwhm_x, fwhm_y, theta_deg) for a 2-D Gaussian - [pixels, pixels, degrees]. + If float/int, this sets the fwhm [pixels] of either Airy disk or + circular Gaussian beam, depending on "mode" parameter. + If 3-element array-like, those are (fwhm_x, fwhm_y, theta_deg) for a 2-D Gaussian + [pixels, pixels, degrees]. In this case, "Airy" mode not supported. if np.array of size > 3, assumes input PSF image - if None, returns original data + if None or 0.0, simply returns original data + mode : str, optional, default "gaussian" + options "gaussian", "airy"; what beam shape to use. Case insensitive. + "airy" only accepts beam of dtype float + or 1-element array-like (i.e., beam must be circular). + this parameter has no effect if beam is a 2-D array. Returns ------- np.array beam-convolved data ''' - if beam is None: + # allow pass through for beam of zero size + if (beam is None): return data - if np.array(beam).size == 1: + + # check inputs + mode = mode.lower() + if mode not in ['gaussian', 'airy']: + raise ValueError(f'mode {mode} not recognized; supported options are "gaussian", "airy"') + if (mode == 'airy') and (np.array(beam).size != 1): + raise ValueError(f'mode "airy" only accepts a single value for the "beam" parameter (distance to first null)') + + if mode == 'airy': + if beam == 0.0: + return data + # convert FWHM to first-null distance + null = 0.5 * (2.44/1.02) * beam + psf = convolution.AiryDisk2DKernel(radius=null) + elif (mode == 'gaussian') and (np.array(beam).size == 1): + if beam == 0.0: + return data fwhm_x = beam fwhm_y = beam theta = 0.0 psf = convolution.Gaussian2DKernel(fwhm_x / 2.35482004503, - fwhm_y / 2.35482004503, - Angle(theta, unit=u.deg)) - elif np.array(beam).size == 3: + fwhm_y / 2.35482004503, + Angle(theta, unit=u.deg)) + elif (mode == 'gaussian') and (np.array(beam).size == 3): (fwhm_x, fwhm_y, theta) = beam psf = convolution.Gaussian2DKernel(fwhm_x / 2.35482004503, - fwhm_y / 2.35482004503, - Angle(theta, unit=u.deg)) + fwhm_y / 2.35482004503, + Angle(theta, unit=u.deg)) else: psf = beam + return convolution.convolve_fft(data, psf) diff --git a/pylanetary/utils/data/Io.yaml b/pylanetary/utils/data/Io.yaml index b1ba6fa..eb6ed11 100644 --- a/pylanetary/utils/data/Io.yaml +++ b/pylanetary/utils/data/Io.yaml @@ -2,8 +2,8 @@ body: name: Io jpl_hor_id: 501 mass: 89.3e21 # kg - req: 1821. # km equatorial radius at 1 bar level - rpol: 1821. # km polar radius + req: 1829.4 # km equatorial radius at 1 bar level + rpol: 1815.7 # km polar radius rvol: 1821. # km volumetric mean radius accel_g: 1.80 # m/s^2 Gravity at 1 bar level num_moons: None From e63d7c44a3256dcdb4c0b71ba5e201610931cec0 Mon Sep 17 00:00:00 2001 From: Edward Molter Date: Thu, 9 Nov 2023 16:24:02 -0800 Subject: [PATCH 2/2] added the ephemeris data for #41 --- pylanetary/navigation/core.py | 27 ++++++++++++++++++- .../navigation/tests/test_nav_remote.py | 4 +-- 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/pylanetary/navigation/core.py b/pylanetary/navigation/core.py index 12c5f46..04bd98f 100755 --- a/pylanetary/navigation/core.py +++ b/pylanetary/navigation/core.py @@ -744,8 +744,31 @@ def write(self, outstem, header={}, flux_unit=''): ''' hdu0 = fits.PrimaryHDU() - hdu0.header = header + # turn ephem into OPAL-like header + # allow to be overwritten by input header dict + header_default = { + 'SIMPLE':True, + 'BITPIX':-32, + 'NAXIS':0, + 'EXTEND':True, + 'NEXTEND':6, + 'TRG_ROT':self.ephem['NPole_ang'], + 'TRG_RA':self.ephem['RA'], + 'TRG_DEC':self.ephem['DEC'], + 'TRG_R_A':self.req, + 'TRG_R_B':self.rpol, + 'TRG_LON':self.ephem['PDObsLon'], + 'TRG_LAT':self.ephem['PDObsLat'], + 'SUN_LAT':self.ephem['PDSunLat'], + 'TRG_PHAS':self.ephem['alpha'], + 'TRG_R':self.ephem['r'], + 'TRG_D':self.ephem['delta'], + } + header_out = {**header_default, **header} #gives priority to input header dict + hdu0.header = fits.Header(header_out) + + # set up metadata for hdul[1-5] hdulist = [hdu0] data_list = [self.data, self.lat_g, @@ -763,6 +786,7 @@ def write(self, outstem, header={}, flux_unit=''): date.format = 'iso' date = date.iso[:10] + # generate hdul[1-5] headers and data for i in range(len(data_list)): data = data_list[i] @@ -779,6 +803,7 @@ def write(self, outstem, header={}, flux_unit=''): hdu = fits.ImageHDU(data=data, header=hdr, name=extnames_list[i]) hdulist.append(hdu) + # write it hdul = fits.HDUList(hdulist) hdul.writeto(outstem, overwrite=True) diff --git a/pylanetary/navigation/tests/test_nav_remote.py b/pylanetary/navigation/tests/test_nav_remote.py index 75f5372..2e9f99e 100644 --- a/pylanetary/navigation/tests/test_nav_remote.py +++ b/pylanetary/navigation/tests/test_nav_remote.py @@ -142,11 +142,11 @@ def test_write(datadir): nav.mu0 = np.cos(np.deg2rad(hdul[5].data)) # write it, then reload it - nav.write(os.path.join(datadir, 'tmp.fits'), header=hdul[0].header, flux_unit='I/F') + nav.write(os.path.join(datadir, 'tmp.fits'), header={}, flux_unit='I/F') hdul_new = fits.open(os.path.join(datadir, 'tmp.fits')) # assert header info propagated - assert hdul_new[0].header == hdul[0].header + assert np.isclose(hdul_new[0].header['TRG_LON'], hdul[0].header['TRG_LON'], rtol=1e-3) # assert data are in correct places for i in range(1,6):