From f84d79bb574336327e37507b4f1311e6f7e70e63 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Thu, 6 Jun 2024 13:10:22 +0100 Subject: [PATCH 1/8] fix fit covariance matrix --- src/adler/science/PhaseCurve.py | 121 ++++++++++++++++++++----- tests/adler/science/test_PhaseCurve.py | 13 ++- 2 files changed, 105 insertions(+), 29 deletions(-) diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index dc53468..c4bfbec 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -1,7 +1,18 @@ from sbpy.photometry import HG, HG1G2, HG12, HG12_Pen16, LinearPhaseFunc import astropy.units as u +import numpy as np from astropy.modeling.fitting import LevMarLSQFitter +# translation between sbpy and adler field names +SBPY_ADLER_DICT = { + "H": "H", + "G": "phase_parameter_1", + "G12": "phase_parameter_1", + "G1": "phase_parameter_1", + "G2": "phase_parameter_2", + "S": "phase_parameter_1", +} + class PhaseCurve: """A class to define the phasecurve model and associated functions. @@ -10,14 +21,14 @@ class PhaseCurve: Attributes ----------- - abs_mag : float + H : float Absolute magnitude, H, of the phasecurve model (often units of mag). - phase_param: float + phase_parameter_1: float The first phase parameter of the phasecurve model, e.g. G from HG (often dimensionless units unless S from LinearPhaseFunc, which has units mag/deg or mag/rad). - phase_param2: float + phase_parameter_2: float The second phase parameter, only used for the 3 parameter HG1G2 phasecurve model. model_name: str @@ -27,30 +38,55 @@ class PhaseCurve: bounds: bool Flag for using the default bounds on the sbpy model (True) or ignoring bounds (False). + H_err : float + Uncertainty in absolute magnitude. + + phase_parameter_1_err : float + Uncertainty in the (first) phase parameter. + + phase_parameter_2_err : float + Uncertainty in the (second, if required) phase parameter. + """ - def __init__(self, abs_mag=18, phase_param=0.2, phase_param2=None, model_name="HG"): - self.abs_mag = abs_mag - self.phase_param = phase_param - self.phase_param2 = phase_param2 + def __init__( + self, + H=18, + phase_parameter_1=0.2, + phase_parameter_2=None, + model_name="HG", + H_err=None, + phase_parameter_1_err=None, + phase_parameter_2_err=None, + ): + self.H = H + self.phase_parameter_1 = phase_parameter_1 + self.phase_parameter_2 = phase_parameter_2 self.model_name = model_name + self.H_err = H_err + self.phase_parameter_1_err = phase_parameter_1_err + self.phase_parameter_2_err = phase_parameter_2_err if model_name == "HG": - self.model_function = HG(H=abs_mag, G=self.phase_param) + self.model_function = HG(H=H, G=self.phase_parameter_1) elif model_name == "HG1G2": - self.model_function = HG1G2(H=abs_mag, G1=self.phase_param, G2=self.phase_param) + self.model_function = HG1G2(H=H, G1=self.phase_parameter_1, G2=self.phase_parameter_1) elif model_name == "HG12": - self.model_function = HG12(H=abs_mag, G12=self.phase_param) + self.model_function = HG12(H=H, G12=self.phase_parameter_1) elif model_name == "HG12_Pen16": - self.model_function = HG12_Pen16(H=abs_mag, G12=self.phase_param) + self.model_function = HG12_Pen16(H=H, G12=self.phase_parameter_1) elif model_name == "LinearPhaseFunc": - self.model_function = LinearPhaseFunc(H=abs_mag, S=self.phase_param) + self.model_function = LinearPhaseFunc(H=H, S=self.phase_parameter_1) else: print("no model selected") + ### set bounds here by default using SetModelBounds? + # for x in model phase parameters: + # self.SetModelBounds(x) + def SetModelBounds(self, param, bound_vals=(None, None)): model_sbpy = self.model_function - param_names = model_sbpy.param_names + # param_names = model_sbpy.param_names x = getattr(model_sbpy, param) setattr(x, "bounds", bound_vals) @@ -74,7 +110,7 @@ def InitModelDict(self, model_dict): Parameters ----------- model_dict : dict - Dictionary containing the PhaseCurve parameters you wish to set, e.g. abs_mag, phase_param + Dictionary containing the PhaseCurve parameters you wish to set, e.g. H, phase_parameter_1 Returns ---------- @@ -109,19 +145,36 @@ def InitModelSbpy(self, model_sbpy): model_name = model_sbpy.__class__.name # get the sbpy model parameters - param_names = model_sbpy.param_names - parameters = [] + param_names = list(model_sbpy.param_names) + # we will also check for any uncertainties we have stored in the sbpy object + param_names_err = ["{}_err".format(x) for x in param_names] + param_names = param_names + param_names_err + + # create a dictionary of phase curve parameters from sbpy in a format accepted by PhaseCurve + parameters = {} for p in param_names: - # try get the quantity (value with units) - x = getattr(model_sbpy, p).quantity - # if there are no units get just the value - if x is None: - x = getattr(model_sbpy, p).value - parameters.append(x) - # print(param_names, parameters) + if p in model_sbpy.__dict__: # check that the parameter is available in the sbpy object + x = getattr(model_sbpy, p) + print(p, x) + # try get the quantity (value with units) + if hasattr(x, "unit"): + print(x.unit) + if (x.unit is None) or ( + x.unit == "" + ): # if there are no units (or weird blank units?) get just the value + x = x.value + else: + x = x.quantity + print(x) + # look up the correct adler parameter name (accounting for additional uncertainty, "_err", parameters) + if p.endswith("_err"): # assumes the uncertainty parameter always ends in "_err" + _p = SBPY_ADLER_DICT[p.split("_err")[0]] + "_err" + parameters[_p] = x + else: + parameters[SBPY_ADLER_DICT[p]] = x # create a PhaseCurve object with the extracted parameters - model = PhaseCurve(*parameters, model_name=model_name) + model = PhaseCurve(**parameters, model_name=model_name) return model @@ -187,6 +240,26 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): else: # unweighted fit model_fit = fitter(self.model_function, phase_angle, reduced_mag) + # Add fitted uncertainties as an additional attribute within the sbpy object + if "param_cov" in fitter.fit_info: + # get the covariance matrix from the fit + covariance = fitter.fit_info["param_cov"] + if covariance is not None: + # get fit uncertainties as square of the diagonal of the covariance + fit_errs = np.sqrt(np.diag(covariance)) + # update only the uncertainties for parameters used in the fit + param_names = np.array(model_fit.param_names) + fit_mask = ~np.array([getattr(model_fit, x).fixed for x in param_names]) + for i, x in enumerate(param_names[fit_mask]): + setattr(model_fit, "{}_err".format(x), fit_errs[i]) + # else: + ### TODO log covariance is None error here + + ### TODO + # else: + # log lack of uncertainties for fitter + # run an MCMC estimate of uncertainty? + ### if overwrite_model: # add an overwrite option? # redo __init__ with the new fitted parameters diff --git a/tests/adler/science/test_PhaseCurve.py b/tests/adler/science/test_PhaseCurve.py index f95b10d..b3c0ceb 100644 --- a/tests/adler/science/test_PhaseCurve.py +++ b/tests/adler/science/test_PhaseCurve.py @@ -10,11 +10,11 @@ def test_PhaseCurve_init(): # test creating a model H = 18.9 G = 0.12 - pc = PhaseCurve(abs_mag=H * u.mag, phase_param=G, model_name="HG") + pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G, model_name="HG") - assert pc.abs_mag.value == 18.9 - assert pc.abs_mag.unit == u.mag - assert pc.phase_param == 0.12 + assert pc.H.value == 18.9 + assert pc.H.unit == u.mag + assert pc.phase_parameter_1 == 0.12 assert pc.model_name == "HG" @@ -27,7 +27,7 @@ def test_PhaseCurve_ReducedMag(): alpha = np.array([0, 10]) * u.deg # linear phase curve model - pc_lin = PhaseCurve(model_name="LinearPhaseFunc", abs_mag=18 * u.mag, phase_param=0.1 * (u.mag / u.deg)) + pc_lin = PhaseCurve(model_name="LinearPhaseFunc", H=18 * u.mag, phase_parameter_1=0.1 * (u.mag / u.deg)) # find the reduced mag red_mag = pc_lin.ReducedMag(alpha) @@ -55,3 +55,6 @@ def test_PhaseCurve_FitModel(): assert pc_fit.H.value == 18.0 assert pc_fit.S.unit == u.mag / u.deg assert pc_fit.S.value == 0.1 + + +### TODO - add more phase curve tests, test bounds, fixed values and fit uncertainties From c1b1c8837fca58af92eb54ada7e17dd3fb341c96 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Thu, 6 Jun 2024 14:35:02 +0100 Subject: [PATCH 2/8] add more tests --- src/adler/science/PhaseCurve.py | 57 +++++++++++++---- tests/adler/science/test_PhaseCurve.py | 87 ++++++++++++++++++++++---- 2 files changed, 121 insertions(+), 23 deletions(-) diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index c4bfbec..99fc16d 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -12,6 +12,14 @@ "G2": "phase_parameter_2", "S": "phase_parameter_1", } +# translation between adler and sbpy fields, depends on model +ADLER_SBPY_DICT = { + "HG": {"phase_parameter_1": "G"}, + "HG1G2": {"phase_parameter_1": "G1", "phase_parameter_2": "G2"}, + "HG12": {"phase_parameter_1": "G12"}, + "HG12_Pen16": {"phase_parameter_1": "G12"}, + "LinearPhaseFunc": {"phase_parameter_1": "S"}, +} class PhaseCurve: @@ -80,15 +88,45 @@ def __init__( else: print("no model selected") - ### set bounds here by default using SetModelBounds? - # for x in model phase parameters: - # self.SetModelBounds(x) - def SetModelBounds(self, param, bound_vals=(None, None)): + """Set the "bounds" attribute of an sbpy model parameter, i.e. the lower and upper constraints for the fitter. + + Parameters + ----------- + param : str + Parameter name bounds to be set fix, if not an sbpy parameter name (e.g. G) the corresponding adler name is looked up (e.g. phase_parameter_1) + + bound_vals : tuple + Set the fitter constraints to (upper, lower) for the param. + + """ model_sbpy = self.model_function - # param_names = model_sbpy.param_names + if param not in model_sbpy.__dict__: + param = ADLER_SBPY_DICT[self.model_name][param] x = getattr(model_sbpy, param) setattr(x, "bounds", bound_vals) + return + + def FixParam(self, param, fix_flag=True): + """Set the "fixed" attribute of an sbpy model parameter. + E.g. use this to fit for absolute magnitude whilst keeping phase parameter fixed. + + Parameters + ----------- + param : str + Parameter name to fix, if not an sbpy parameter name (e.g. G) the corresponding adler name is looked up (e.g. phase_parameter_1) + + fix_flag : bool + Set True to keep the param fixed when fitting + + """ + + model_sbpy = self.model_function + if param not in model_sbpy.__dict__: + param = ADLER_SBPY_DICT[self.model_name][param] + x = getattr(model_sbpy, param) + setattr(x, "fixed", fix_flag) + return def ReturnModelDict(self): """Return the values for the PhaseCurve class as a dict @@ -127,7 +165,7 @@ def InitModelDict(self, model_dict): def InitModelSbpy(self, model_sbpy): """Set up a new PhaseCurve model object from an existing sbpy model - ### or create dict from sbpy model and then use InitModelDict? + ### TODO or create dict from sbpy model and then use InitModelDict? Parameters ----------- @@ -155,17 +193,14 @@ def InitModelSbpy(self, model_sbpy): for p in param_names: if p in model_sbpy.__dict__: # check that the parameter is available in the sbpy object x = getattr(model_sbpy, p) - print(p, x) # try get the quantity (value with units) if hasattr(x, "unit"): - print(x.unit) if (x.unit is None) or ( x.unit == "" ): # if there are no units (or weird blank units?) get just the value x = x.value else: x = x.quantity - print(x) # look up the correct adler parameter name (accounting for additional uncertainty, "_err", parameters) if p.endswith("_err"): # assumes the uncertainty parameter always ends in "_err" _p = SBPY_ADLER_DICT[p.split("_err")[0]] + "_err" @@ -220,7 +255,8 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): Uncertainty on the reduced magnitude, used to weight the fit. fitter : object - Select a fitting function from astropy.modeling.fitting, defaults to astropy.modeling.fitting.LevMarLSQFitter + Select a fitting function from astropy.modeling.fitting, defaults to astropy.modeling.fitting.LevMarLSQFitter. + N.B. that LevMarLSQFitter cannot handle inequality constraints for the HG1G2 model. Returns ---------- @@ -233,7 +269,6 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): # use the LevMarLSQFitter by default if fitter is None: fitter = LevMarLSQFitter() - # print(fitter) if mag_err is not None: # fit weighted by photometric uncertainty model_fit = fitter(self.model_function, phase_angle, reduced_mag, weights=1.0 / mag_err) diff --git a/tests/adler/science/test_PhaseCurve.py b/tests/adler/science/test_PhaseCurve.py index b3c0ceb..4cbf76b 100644 --- a/tests/adler/science/test_PhaseCurve.py +++ b/tests/adler/science/test_PhaseCurve.py @@ -1,11 +1,12 @@ +from adler.science.PhaseCurve import PhaseCurve from numpy.testing import assert_array_equal import pytest +import numpy as np +import astropy.units as u def test_PhaseCurve_init(): - import numpy as np - import astropy.units as u - from adler.science.PhaseCurve import PhaseCurve + """Test intialising the model.""" # test creating a model H = 18.9 @@ -19,9 +20,7 @@ def test_PhaseCurve_init(): def test_PhaseCurve_ReducedMag(): - import numpy as np - import astropy.units as u - from adler.science.PhaseCurve import PhaseCurve + """Test calculating the reduced magnitude from a PhaseCurve object.""" # define the phase angles alpha = np.array([0, 10]) * u.deg @@ -36,10 +35,8 @@ def test_PhaseCurve_ReducedMag(): assert_array_equal(red_mag.value, np.array([18.0, 19.0])) -def test_PhaseCurve_FitModel(): - import numpy as np - import astropy.units as u - from adler.science.PhaseCurve import PhaseCurve +def test_PhaseCurve_FitModel_linear(): + """Test fitting a linear phase function to two data points.""" # define the observations alpha = np.array([0, 10]) * u.deg @@ -48,13 +45,79 @@ def test_PhaseCurve_FitModel(): # empty linear phase curve model pc_lin = PhaseCurve(model_name="LinearPhaseFunc") - # fit the model to the data + # fit the model to the data - this returns an sbpy object pc_fit = pc_lin.FitModel(alpha, red_mag) assert pc_fit.H.unit == u.mag assert pc_fit.H.value == 18.0 assert pc_fit.S.unit == u.mag / u.deg assert pc_fit.S.value == 0.1 + assert ( + hasattr(pc_fit, "H_err") == False + ) # with only two data points the covariance matrix is None and no uncertainties are stored -### TODO - add more phase curve tests, test bounds, fixed values and fit uncertainties +def test_PhaseCurve_FitModel_HG(): + """Test fitting a HG model to generated data.""" + + # generate some model data + pc1 = PhaseCurve(H=18.0 * u.mag, phase_parameter_1=0.15, model_name="HG") + alpha = np.linspace(0, 30) * u.deg + red_mag = pc1.ReducedMag(alpha) + + # fit the same phase curve model to the data + pc_fit = pc1.FitModel(alpha, red_mag) + # convert from sbpy to adler PhaseCurve object + pc2 = pc1.InitModelSbpy(pc_fit) + + # the new fitted model should have the same parameters as the input model + assert pc2.H == pc1.H + assert pc2.phase_parameter_1 == pc1.phase_parameter_1 + assert pc2.phase_parameter_2 is None + assert pc1.phase_parameter_1_err is None # the first model had no uncertainties + assert pc2.phase_parameter_1_err is not None # the fitted model has some uncertainties + + +def test_PhaseCurve_FitModel_HG_fixed(): + """Test fitting a just H whilst keeping G fixed.""" + + # generate some model data + pc1 = PhaseCurve(H=18.0 * u.mag, phase_parameter_1=0.15, model_name="HG") + alpha = np.linspace(0, 30) * u.deg + red_mag = pc1.ReducedMag(alpha) + + # fix phase_parameter_1 + pc1.FixParam("phase_parameter_1") + + # fit the same phase curve model to the data + pc_fit = pc1.FitModel(alpha, red_mag) + # convert from sbpy to adler PhaseCurve object + pc2 = pc1.InitModelSbpy(pc_fit) + + # the new fitted model should have the same parameters as the input model, but G is fixed + assert pc_fit._fixed["G"] is True + assert pc2.H == pc1.H + assert pc2.phase_parameter_1 == pc1.phase_parameter_1 + assert pc2.phase_parameter_1_err is None # the fitted model has no uncertainties when param is fixed + + +def test_PhaseCurve_FitModel_HG_bounds(): + """Test fitting a just H whilst keeping G fixed.""" + + # generate some model data + pc1 = PhaseCurve(H=18.0 * u.mag, phase_parameter_1=0.15, model_name="HG") + alpha = np.linspace(0, 30) * u.deg + red_mag = pc1.ReducedMag(alpha) + + # set bounds on phase parameter + pc1.SetModelBounds("phase_parameter_1", (0.0, 0.1)) + + # fit the same phase curve model to the data + pc_fit = pc1.FitModel(alpha, red_mag) + # convert from sbpy to adler PhaseCurve object + pc2 = pc1.InitModelSbpy(pc_fit) + + # the new fitted model should have the same parameters as the input model, but G is fixed + assert pc_fit.G.bounds == (0.0, 0.1) + assert pc2.phase_parameter_1 == 0.1 + assert pc2.phase_parameter_1_err is not None From 2de2e322c8a052323c53e1cd936681d53fd4ef53 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Thu, 6 Jun 2024 15:41:55 +0100 Subject: [PATCH 3/8] no units test --- src/adler/science/PhaseCurve.py | 3 ++- tests/adler/science/test_PhaseCurve.py | 22 ++++++++++++++++++++++ 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index 99fc16d..a40d1cb 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -256,7 +256,7 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): fitter : object Select a fitting function from astropy.modeling.fitting, defaults to astropy.modeling.fitting.LevMarLSQFitter. - N.B. that LevMarLSQFitter cannot handle inequality constraints for the HG1G2 model. + N.B. that LevMarLSQFitter cannot handle inequality constraints for the HG1G2 model, use something like SLSQPLSQFitter from astropy.modeling.fitting (does not return covariance matrix!). Returns ---------- @@ -297,5 +297,6 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): ### if overwrite_model: # add an overwrite option? # redo __init__ with the new fitted parameters + # this would then return an adler PhaseCurve object rather than an sbpy object return model_fit diff --git a/tests/adler/science/test_PhaseCurve.py b/tests/adler/science/test_PhaseCurve.py index 4cbf76b..77e0385 100644 --- a/tests/adler/science/test_PhaseCurve.py +++ b/tests/adler/science/test_PhaseCurve.py @@ -78,6 +78,28 @@ def test_PhaseCurve_FitModel_HG(): assert pc2.phase_parameter_1_err is not None # the fitted model has some uncertainties +def test_PhaseCurve_FitModel_HG_no_units(): + """Test fitting a HG model to generated data, but without units. + If units are not provided, the phase angles must be in radians!""" + + # generate some model data + pc1 = PhaseCurve(H=18.0, phase_parameter_1=0.15, model_name="HG") + alpha = np.radians(np.linspace(0, 30)) + red_mag = pc1.ReducedMag(alpha) + + # fit the same phase curve model to the data + pc_fit = pc1.FitModel(alpha, red_mag) + # convert from sbpy to adler PhaseCurve object + pc2 = pc1.InitModelSbpy(pc_fit) + + # the new fitted model should have the same parameters as the input model + assert pc2.H == pc1.H + assert pc2.phase_parameter_1 == pc1.phase_parameter_1 + assert pc2.phase_parameter_2 is None + assert pc1.phase_parameter_1_err is None # the first model had no uncertainties + assert pc2.phase_parameter_1_err is not None # the fitted model has some uncertainties + + def test_PhaseCurve_FitModel_HG_fixed(): """Test fitting a just H whilst keeping G fixed.""" From a224c194165d3e75709f4c167844b8fa840055e4 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Thu, 6 Jun 2024 15:52:17 +0100 Subject: [PATCH 4/8] update PhaseCurve parameter names --- src/adler/adler.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/adler/adler.py b/src/adler/adler.py index 2987cc6..4cbca5d 100644 --- a/src/adler/adler.py +++ b/src/adler/adler.py @@ -100,8 +100,8 @@ def runAdler(cli_args): # initial simple phase curve filter model with fixed G12 pc = PhaseCurve( - abs_mag=sso.H * u.mag, - phase_param=0.62, + H=sso.H * u.mag, + phase_parameter_1=0.62, model_name="HG12_Pen16", ) @@ -168,15 +168,15 @@ def runAdler(cli_args): # if sum(obs_mask) < N_pc_fit: # # use an assumed value of G12 until more data is available - # pc_fit = PhaseCurve(abs_mag=sso.H * u.mag, phase_param=0.62, model_name="HG12_Pen16") + # pc_fit = PhaseCurve(H=sso.H * u.mag, phase_parameter_1=0.62, model_name="HG12_Pen16") # else: # # do a simple HG12_Pen16 fit to the past data # pc_fit = pc.FitModel(alpha[obs_mask], red_mag[obs_mask], mag_err[obs_mask]) # pc_fit = pc.InitModelSbpy(pc_fit) # print(pc_fit) - # print(pc_fit.abs_mag) - # print(pc_fit.phase_param) + # print(pc_fit.H) + # print(pc_fit.phase_parameter_1) # # now check if the new observations are outlying # alpha_new = alpha[~obs_mask] From 6f067fc514b7588c3184c971ad3b41ead55207a7 Mon Sep 17 00:00:00 2001 From: Steph Merritt <97111051+astronomerritt@users.noreply.github.com> Date: Tue, 18 Jun 2024 11:35:33 +0100 Subject: [PATCH 5/8] Adding JSON constructor methods. (#144) --- src/adler/dataclasses/AdlerPlanetoid.py | 28 ++++++++++++++++++++ src/adler/dataclasses/MPCORB.py | 11 +++++++- src/adler/dataclasses/Observations.py | 15 ++++++++++- src/adler/dataclasses/SSObject.py | 23 +++++++++++++++- src/adler/dataclasses/dataclass_utilities.py | 23 +++++++++++++++- 5 files changed, 96 insertions(+), 4 deletions(-) diff --git a/src/adler/dataclasses/AdlerPlanetoid.py b/src/adler/dataclasses/AdlerPlanetoid.py index 18748bd..d6ceb40 100644 --- a/src/adler/dataclasses/AdlerPlanetoid.py +++ b/src/adler/dataclasses/AdlerPlanetoid.py @@ -1,6 +1,8 @@ from lsst.rsp import get_tap_service import pandas as pd +import numpy as np import logging +import json from adler.dataclasses.Observations import Observations from adler.dataclasses.MPCORB import MPCORB @@ -114,6 +116,32 @@ def construct_from_SQL( return cls(ssObjectId, filter_list, date_range, observations_by_filter, mpcorb, ssobject, adler_data) + @classmethod + def construct_from_JSON(cls, json_filename): + with open(json_filename) as f: + json_dict = json.load(f) + + observations_dict = {**json_dict["SSSource"], **json_dict["DiaSource"]} + + filter_list = [observations_dict["band"]] + + MPCORB_dict = json_dict["MPCORB"] + SSObject_dict = json_dict["SSObject"] + + ssObjectId = observations_dict["ssObjectId"] + + observations_by_filter = [ + Observations.construct_from_dictionary(ssObjectId, filter_list[0], observations_dict) + ] + mpcorb = MPCORB.construct_from_dictionary(ssObjectId, MPCORB_dict) + ssobject = SSObject.construct_from_dictionary(ssObjectId, filter_list, SSObject_dict) + + adler_data = AdlerData(ssObjectId, filter_list) + + return cls( + ssObjectId, filter_list, [np.nan, np.nan], observations_by_filter, mpcorb, ssobject, adler_data + ) + @classmethod def construct_from_RSP( cls, ssObjectId, filter_list=["u", "g", "r", "i", "z", "y"], date_range=[60000.0, 67300.0] diff --git a/src/adler/dataclasses/MPCORB.py b/src/adler/dataclasses/MPCORB.py index a17bd18..289084d 100644 --- a/src/adler/dataclasses/MPCORB.py +++ b/src/adler/dataclasses/MPCORB.py @@ -1,6 +1,6 @@ from dataclasses import dataclass -from adler.dataclasses.dataclass_utilities import get_from_table +from adler.dataclasses.dataclass_utilities import get_from_table, get_from_dictionary MPCORB_KEYS = { "mpcDesignation": str, @@ -109,3 +109,12 @@ def construct_from_data_table(cls, ssObjectId, data_table): mpcorb_dict[mpcorb_key] = get_from_table(data_table, mpcorb_key, mpcorb_type, "MPCORB") return cls(**mpcorb_dict) + + @classmethod + def construct_from_dictionary(cls, ssObjectId, data_dict): + mpcorb_dict = {"ssObjectId": ssObjectId} + + for mpcorb_key, mpcorb_type in MPCORB_KEYS.items(): + mpcorb_dict[mpcorb_key] = get_from_dictionary(data_dict, mpcorb_key, mpcorb_type, "MPCORB") + + return cls(**mpcorb_dict) diff --git a/src/adler/dataclasses/Observations.py b/src/adler/dataclasses/Observations.py index 64b3900..e99a758 100644 --- a/src/adler/dataclasses/Observations.py +++ b/src/adler/dataclasses/Observations.py @@ -1,7 +1,7 @@ from dataclasses import dataclass, field import numpy as np -from adler.dataclasses.dataclass_utilities import get_from_table +from adler.dataclasses.dataclass_utilities import get_from_table, get_from_dictionary OBSERVATIONS_KEYS = { "mag": np.ndarray, @@ -107,6 +107,19 @@ def construct_from_data_table(cls, ssObjectId, filter_name, data_table): return cls(**obs_dict) + @classmethod + def construct_from_dictionary(cls, ssObjectId, filter_name, data_dict): + obs_dict = {"ssObjectId": ssObjectId, "filter_name": filter_name, "num_obs": 1} + + for obs_key, obs_type in OBSERVATIONS_KEYS.items(): + obs_dict[obs_key] = get_from_dictionary(data_dict, obs_key, obs_type, "SSSource/DIASource") + + obs_dict["reduced_mag"] = cls.calculate_reduced_mag( + cls, obs_dict["mag"], obs_dict["topocentricDist"], obs_dict["heliocentricDist"] + ) + + return cls(**obs_dict) + def calculate_reduced_mag(self, mag, topocentric_dist, heliocentric_dist): """ Calculates the reduced magnitude column. diff --git a/src/adler/dataclasses/SSObject.py b/src/adler/dataclasses/SSObject.py index 9ec0443..eedc724 100644 --- a/src/adler/dataclasses/SSObject.py +++ b/src/adler/dataclasses/SSObject.py @@ -1,7 +1,7 @@ from dataclasses import dataclass, field import numpy as np -from adler.dataclasses.dataclass_utilities import get_from_table +from adler.dataclasses.dataclass_utilities import get_from_table, get_from_dictionary SSO_KEYS = { "discoverySubmissionDate": float, @@ -86,6 +86,27 @@ def construct_from_data_table(cls, ssObjectId, filter_list, data_table): return cls(**sso_dict) + @classmethod + def construct_from_dictionary(cls, ssObjectId, filter_list, data_dict): + sso_dict = {"ssObjectId": ssObjectId, "filter_list": filter_list, "filter_dependent_values": []} + + for sso_key, sso_type in SSO_KEYS.items(): + sso_dict[sso_key] = get_from_dictionary(data_dict, sso_key, sso_type, "SSObject") + + for i, filter_name in enumerate(filter_list): + filter_dept_object = FilterDependentSSO( + filter_name=filter_name, + H=get_from_dictionary(data_dict, filter_name + "_H", float, "SSObject"), + G12=get_from_dictionary(data_dict, filter_name + "_G12", float, "SSObject"), + Herr=get_from_dictionary(data_dict, filter_name + "_HErr", float, "SSObject"), + G12err=get_from_dictionary(data_dict, filter_name + "_G12Err", float, "SSObject"), + nData=get_from_dictionary(data_dict, filter_name + "_Ndata", float, "SSObject"), + ) + + sso_dict["filter_dependent_values"].append(filter_dept_object) + + return cls(**sso_dict) + @dataclass class FilterDependentSSO: diff --git a/src/adler/dataclasses/dataclass_utilities.py b/src/adler/dataclasses/dataclass_utilities.py index 40d5818..d82c031 100644 --- a/src/adler/dataclasses/dataclass_utilities.py +++ b/src/adler/dataclasses/dataclass_utilities.py @@ -86,7 +86,7 @@ def get_from_table(data_table, column_name, data_type, table_name="default"): elif data_type == int: data_val = int(data_table[column_name][0]) elif data_type == np.ndarray: - data_val = np.array(data_table[column_name]) + data_val = np.array(data_table[column_name], ndmin=1) else: logger.error( "TypeError: Type for argument data_type not recognised for column {} in table {}: must be str, float, int or np.ndarray.".format( @@ -108,6 +108,27 @@ def get_from_table(data_table, column_name, data_type, table_name="default"): return data_val +def get_from_dictionary(data_dict, key_name, data_type, table_name="default"): + try: + if data_type == str: + data_val = str(data_dict[key_name]) + elif data_type == float: + data_val = float(data_dict[key_name]) + elif data_type == int: + data_val = int(data_dict[key_name]) + elif data_type == np.ndarray: + data_val = np.array(data_dict[key_name], ndmin=1) + else: + print("type not recognised") + + except ValueError: + print("error message") + + data_val = check_value_populated(data_val, data_type, key_name, "JSON") + + return data_val + + def check_value_populated(data_val, data_type, column_name, table_name): """Checks to see if data_val populated properly and prints a helpful warning if it didn't. Usually this will trigger because the RSP hasn't populated that field for this particular object. From 1ccc7109732f26deb59c81d04e16a8e8769e7ffc Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 18 Jun 2024 19:04:08 +0000 Subject: [PATCH 6/8] ephem test nb --- notebooks/simulated_asteroid_ephem_rsp.ipynb | 657 +++++++++++++++++++ 1 file changed, 657 insertions(+) create mode 100644 notebooks/simulated_asteroid_ephem_rsp.ipynb diff --git a/notebooks/simulated_asteroid_ephem_rsp.ipynb b/notebooks/simulated_asteroid_ephem_rsp.ipynb new file mode 100644 index 0000000..7c4b0df --- /dev/null +++ b/notebooks/simulated_asteroid_ephem_rsp.ipynb @@ -0,0 +1,657 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Testing DP0.3 ephemerides and coordinates\n", + "By Jamie Robinson" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook retrieves the observations and orbital parameters of an object in DP0.3. We compare the ephemerides and coordinates from DP0.3 to the output from propagating the orbit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:15.026373Z", + "iopub.status.busy": "2024-06-18T19:01:15.026079Z", + "iopub.status.idle": "2024-06-18T19:01:17.359920Z", + "shell.execute_reply": "2024-06-18T19:01:17.359236Z", + "shell.execute_reply.started": "2024-06-18T19:01:15.026354Z" + } + }, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "from astropy import units as u\n", + "from astropy.coordinates import SkyCoord, GCRS\n", + "from sbpy.data import Orbit\n", + "from astropy.time import Time\n", + "from astropy.coordinates import solar_system_ephemeris\n", + "from sbpy.photometry import HG\n", + "from astropy.table import QTable\n", + "\n", + "from lsst.rsp import get_tap_service" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:17.361390Z", + "iopub.status.busy": "2024-06-18T19:01:17.360891Z", + "iopub.status.idle": "2024-06-18T19:01:17.438804Z", + "shell.execute_reply": "2024-06-18T19:01:17.437930Z", + "shell.execute_reply.started": "2024-06-18T19:01:17.361368Z" + } + }, + "outputs": [], + "source": [ + "service = get_tap_service(\"ssotap\")\n", + "assert service is not None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:17.439847Z", + "iopub.status.busy": "2024-06-18T19:01:17.439631Z", + "iopub.status.idle": "2024-06-18T19:01:17.442630Z", + "shell.execute_reply": "2024-06-18T19:01:17.442049Z", + "shell.execute_reply.started": "2024-06-18T19:01:17.439830Z" + } + }, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "ssoid = \"6098332225018\" # good test object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:17.443664Z", + "iopub.status.busy": "2024-06-18T19:01:17.443428Z", + "iopub.status.idle": "2024-06-18T19:01:18.508689Z", + "shell.execute_reply": "2024-06-18T19:01:18.508031Z", + "shell.execute_reply.started": "2024-06-18T19:01:17.443646Z" + } + }, + "outputs": [], + "source": [ + "query = \"\"\"\n", + "SELECT\n", + " *\n", + "FROM\n", + " dp03_catalogs_10yr.DiaSource as dia\n", + "INNER JOIN\n", + " dp03_catalogs_10yr.SSSource as sss\n", + "ON\n", + " dia.diaSourceId = sss.diaSourceId\n", + "WHERE\n", + " dia.ssObjectId={}\n", + "ORDER by dia.ssObjectId\n", + "\"\"\".format(\n", + " ssoid\n", + ")\n", + "\n", + "df = service.search(query).to_table().to_pandas()\n", + "df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.509807Z", + "iopub.status.busy": "2024-06-18T19:01:18.509560Z", + "iopub.status.idle": "2024-06-18T19:01:18.661531Z", + "shell.execute_reply": "2024-06-18T19:01:18.660953Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.509790Z" + } + }, + "outputs": [], + "source": [ + "results = service.search(\"SELECT * FROM dp03_catalogs_10yr.MPCORB \" \"WHERE ssObjectId={}\".format(ssoid))\n", + "df_orb = results.to_table().to_pandas()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.663655Z", + "iopub.status.busy": "2024-06-18T19:01:18.663419Z", + "iopub.status.idle": "2024-06-18T19:01:18.683552Z", + "shell.execute_reply": "2024-06-18T19:01:18.682809Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.663638Z" + } + }, + "outputs": [], + "source": [ + "df_orb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.685016Z", + "iopub.status.busy": "2024-06-18T19:01:18.684704Z", + "iopub.status.idle": "2024-06-18T19:01:18.699123Z", + "shell.execute_reply": "2024-06-18T19:01:18.698473Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.684988Z" + } + }, + "outputs": [], + "source": [ + "# put obs in time order\n", + "df = df.sort_values(\"midPointMjdTai\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Resample the sparse observations by propagating orbit with sbpy & oorb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.700389Z", + "iopub.status.busy": "2024-06-18T19:01:18.700015Z", + "iopub.status.idle": "2024-06-18T19:01:18.715818Z", + "shell.execute_reply": "2024-06-18T19:01:18.715179Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.700367Z" + } + }, + "outputs": [], + "source": [ + "# check orbit\n", + "e = df_orb[\"e\"]\n", + "incl = df_orb[\"incl\"]\n", + "q = df_orb[\"q\"]\n", + "a = q / (1.0 - e)\n", + "Q = a * (1.0 + e)\n", + "print(a, e, incl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.716909Z", + "iopub.status.busy": "2024-06-18T19:01:18.716682Z", + "iopub.status.idle": "2024-06-18T19:01:18.743743Z", + "shell.execute_reply": "2024-06-18T19:01:18.743021Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.716892Z" + } + }, + "outputs": [], + "source": [ + "df_orb[\"a\"] = a\n", + "df_orb[\"Q\"] = Q\n", + "df_orb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.745120Z", + "iopub.status.busy": "2024-06-18T19:01:18.744684Z", + "iopub.status.idle": "2024-06-18T19:01:18.761040Z", + "shell.execute_reply": "2024-06-18T19:01:18.760412Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.745100Z" + } + }, + "outputs": [], + "source": [ + "# Calculate some extra orbital elements\n", + "df_orb[\"P\"] = df_orb[\"a\"] ** (3.0 / 2.0) # orbital period in years\n", + "df_orb[\"n\"] = 360.0 / (df_orb[\"P\"] * 365.25) # mean motion in deg/day\n", + "df_orb[\"M\"] = (\n", + " df_orb[\"n\"] * (df_orb[\"epoch\"] - df_orb[\"tperi\"])\n", + ") % 360 # angles must be in correct range otherwise sbpy/pyoorb freak out\n", + "df_orb[[\"a\", \"P\", \"n\", \"M\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.762022Z", + "iopub.status.busy": "2024-06-18T19:01:18.761789Z", + "iopub.status.idle": "2024-06-18T19:01:18.775973Z", + "shell.execute_reply": "2024-06-18T19:01:18.775117Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.762000Z" + } + }, + "outputs": [], + "source": [ + "# rename columns for consistency\n", + "df_orb = df_orb.rename(columns={\"node\": \"Omega\", \"peri\": \"w\"})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.777363Z", + "iopub.status.busy": "2024-06-18T19:01:18.776889Z", + "iopub.status.idle": "2024-06-18T19:01:18.796868Z", + "shell.execute_reply": "2024-06-18T19:01:18.796165Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.777340Z" + } + }, + "outputs": [], + "source": [ + "df_orb[[\"a\", \"e\", \"incl\", \"Omega\", \"w\", \"M\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.798419Z", + "iopub.status.busy": "2024-06-18T19:01:18.797782Z", + "iopub.status.idle": "2024-06-18T19:01:18.816304Z", + "shell.execute_reply": "2024-06-18T19:01:18.815461Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.798397Z" + } + }, + "outputs": [], + "source": [ + "# create an sbpy oorb object from dataframe via QTable\n", + "tab = QTable.from_pandas(\n", + " df_orb[[\"a\", \"e\", \"incl\", \"Omega\", \"w\", \"M\"]],\n", + " units={\"a\": u.au, \"incl\": u.deg, \"Omega\": u.deg, \"w\": u.deg, \"M\": u.deg},\n", + ")\n", + "orbit = Orbit.from_table(tab)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.817772Z", + "iopub.status.busy": "2024-06-18T19:01:18.817316Z", + "iopub.status.idle": "2024-06-18T19:01:18.833560Z", + "shell.execute_reply": "2024-06-18T19:01:18.832844Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.817749Z" + } + }, + "outputs": [], + "source": [ + "# oorb requires certain extra fields\n", + "orbit[\"epoch\"] = Time(Time(df_orb[\"epoch\"], format=\"mjd\").jd, format=\"jd\")\n", + "orbit[\"targetname\"] = np.array(df_orb[\"ssObjectId\"]).astype(str)\n", + "orbit[\"H\"] = df_orb[\"mpcH\"] * u.mag\n", + "orbit[\"G\"] = df_orb[\"mpcG\"] * u.dimensionless_unscaled" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.834791Z", + "iopub.status.busy": "2024-06-18T19:01:18.834470Z", + "iopub.status.idle": "2024-06-18T19:01:18.848589Z", + "shell.execute_reply": "2024-06-18T19:01:18.847950Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.834772Z" + } + }, + "outputs": [], + "source": [ + "orbit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.849539Z", + "iopub.status.busy": "2024-06-18T19:01:18.849320Z", + "iopub.status.idle": "2024-06-18T19:01:18.865983Z", + "shell.execute_reply": "2024-06-18T19:01:18.865340Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.849513Z" + } + }, + "outputs": [], + "source": [ + "# define a set of JD times to propagate the orbital elements to\n", + "N = 1000\n", + "times = Time(\n", + " Time(np.linspace(np.amin(df[\"midPointMjdTai\"]), np.amax(df[\"midPointMjdTai\"]), N), format=\"mjd\").jd,\n", + " format=\"jd\",\n", + ")\n", + "times[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.867204Z", + "iopub.status.busy": "2024-06-18T19:01:18.866796Z", + "iopub.status.idle": "2024-06-18T19:01:18.880152Z", + "shell.execute_reply": "2024-06-18T19:01:18.879521Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.867185Z" + } + }, + "outputs": [], + "source": [ + "# create an empty dataframe to hold resampled observations\n", + "df_dense = pd.DataFrame()\n", + "df_dense[\"midPointMjdTai\"] = times.mjd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.881123Z", + "iopub.status.busy": "2024-06-18T19:01:18.880917Z", + "iopub.status.idle": "2024-06-18T19:01:51.982325Z", + "shell.execute_reply": "2024-06-18T19:01:51.981730Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.881108Z" + } + }, + "outputs": [], + "source": [ + "# propagate the orbit forward in time.\n", + "# probably a better way to do this but I can't get oo_propagate to work with a Time list right now\n", + "# see: https://github.com/NASA-Planetary-Science/sbpy/issues/341\n", + "\n", + "df_pos = pd.DataFrame() # empty dataframe to hold cartesian coordinates\n", + "\n", + "for i in range(len(times)):\n", + " print(i)\n", + " prop_elem = orbit.oo_propagate(times[i]) # propagate the orbit to the selected time step\n", + " del prop_elem.table[\n", + " \"orbtype\"\n", + " ] # orbtype is added as int, sbpy freaks out so delete the orbtype and then _to_oo works it out\n", + " print(\"propagate\")\n", + " statevec = prop_elem.oo_transform(\"CART\") # transform from orbital elements to cartesian\n", + " print(\"transform\")\n", + "\n", + " # append new cartesian coordinates to the dataframe\n", + " _df_statevec = statevec.table.to_pandas()\n", + " df_pos = pd.concat((df_pos, _df_statevec))\n", + "\n", + "df_pos.reset_index(drop=True, inplace=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:51.983360Z", + "iopub.status.busy": "2024-06-18T19:01:51.983151Z", + "iopub.status.idle": "2024-06-18T19:01:52.000157Z", + "shell.execute_reply": "2024-06-18T19:01:51.999473Z", + "shell.execute_reply.started": "2024-06-18T19:01:51.983345Z" + } + }, + "outputs": [], + "source": [ + "df_pos" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Use astropy coordinates to transform between all the coordinate systems" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.001741Z", + "iopub.status.busy": "2024-06-18T19:01:52.001163Z", + "iopub.status.idle": "2024-06-18T19:01:52.010357Z", + "shell.execute_reply": "2024-06-18T19:01:52.009551Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.001719Z" + } + }, + "outputs": [], + "source": [ + "# define heliocentric cartesian coordinates\n", + "c_xyz_hel = SkyCoord(\n", + " x=np.array(df_pos[\"x\"]),\n", + " y=np.array(df_pos[\"y\"]),\n", + " z=np.array(df_pos[\"z\"]),\n", + " unit=\"AU\",\n", + " representation_type=\"cartesian\",\n", + " frame=\"heliocentrictrueecliptic\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.013233Z", + "iopub.status.busy": "2024-06-18T19:01:52.012889Z", + "iopub.status.idle": "2024-06-18T19:01:52.025365Z", + "shell.execute_reply": "2024-06-18T19:01:52.024660Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.013210Z" + } + }, + "outputs": [], + "source": [ + "# transform to heliocentric ecliptic coords\n", + "c_ecl_hel = c_xyz_hel.copy()\n", + "c_ecl_hel.representation_type = \"spherical\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.026475Z", + "iopub.status.busy": "2024-06-18T19:01:52.026241Z", + "iopub.status.idle": "2024-06-18T19:01:52.095939Z", + "shell.execute_reply": "2024-06-18T19:01:52.095204Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.026457Z" + } + }, + "outputs": [], + "source": [ + "# transform to geocentric equatorial coords (times required to calculate Earth position)\n", + "with solar_system_ephemeris.set(\"jpl\"):\n", + " c_eq_geo = c_xyz_hel.transform_to(GCRS(obstime=times))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.096956Z", + "iopub.status.busy": "2024-06-18T19:01:52.096748Z", + "iopub.status.idle": "2024-06-18T19:01:52.100730Z", + "shell.execute_reply": "2024-06-18T19:01:52.099911Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.096940Z" + } + }, + "outputs": [], + "source": [ + "# transform to geocentric cartesian coords\n", + "c_xyz_geo = c_eq_geo.copy()\n", + "c_xyz_geo.representation_type = \"cartesian\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.101667Z", + "iopub.status.busy": "2024-06-18T19:01:52.101450Z", + "iopub.status.idle": "2024-06-18T19:01:52.118900Z", + "shell.execute_reply": "2024-06-18T19:01:52.118257Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.101652Z" + } + }, + "outputs": [], + "source": [ + "# transform from geo equatorial (ra, dec) to geo ecliptic (lon, lat)\n", + "c_ecl_geo = c_eq_geo.transform_to(\"geocentrictrueecliptic\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.120101Z", + "iopub.status.busy": "2024-06-18T19:01:52.119878Z", + "iopub.status.idle": "2024-06-18T19:01:52.364967Z", + "shell.execute_reply": "2024-06-18T19:01:52.364185Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.120086Z" + } + }, + "outputs": [], + "source": [ + "# plot the propagated cartesian positions against the database values\n", + "\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot1 = \"heliocentricX\"\n", + "y_plot2 = \"heliocentricY\"\n", + "y_plot3 = \"heliocentricZ\"\n", + "df_plot = df\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot1], label=y_plot1)\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot2], label=y_plot2)\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot3], label=y_plot3)\n", + "\n", + "ax1.plot(times.mjd, c_xyz_hel.x)\n", + "ax1.plot(times.mjd, c_xyz_hel.y)\n", + "ax1.plot(times.mjd, c_xyz_hel.z)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(\"distance\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are some deviations of x and y positions, probably due to slightly different reference frames and methods of propagating orbits.\n", + "\n", + "There is something wrong with database z positions, will be fixed soon!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.366083Z", + "iopub.status.busy": "2024-06-18T19:01:52.365860Z", + "iopub.status.idle": "2024-06-18T19:01:52.558412Z", + "shell.execute_reply": "2024-06-18T19:01:52.557808Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.366066Z" + } + }, + "outputs": [], + "source": [ + "# the ecliptic coordinates look good!\n", + "\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot1 = \"eclipticLambda\"\n", + "y_plot2 = \"eclipticBeta\"\n", + "df_plot = df\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot1], label=y_plot1)\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot2], label=y_plot2)\n", + "\n", + "ax1.plot(times.mjd, c_ecl_geo.lon.degree)\n", + "ax1.plot(times.mjd, c_ecl_geo.lat.degree)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(\"angle\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "lsst-sbpy", + "language": "python", + "name": "lsst-sbpy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 922a5404a332e001b6074ce5720aa14940333abb Mon Sep 17 00:00:00 2001 From: jrob93 Date: Fri, 5 Jul 2024 15:41:05 +0000 Subject: [PATCH 7/8] nbs comparing dp0.3 eph to horizons/astropy eph --- .../simulated_asteroid_ephem_rsp.ipynb | 238 +++-- ...imulated_asteroid_ephem_rsp_horizons.ipynb | 910 ++++++++++++++++++ 2 files changed, 1018 insertions(+), 130 deletions(-) rename notebooks/{ => sim_ephem}/simulated_asteroid_ephem_rsp.ipynb (61%) create mode 100644 notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb diff --git a/notebooks/simulated_asteroid_ephem_rsp.ipynb b/notebooks/sim_ephem/simulated_asteroid_ephem_rsp.ipynb similarity index 61% rename from notebooks/simulated_asteroid_ephem_rsp.ipynb rename to notebooks/sim_ephem/simulated_asteroid_ephem_rsp.ipynb index 7c4b0df..b890d7f 100644 --- a/notebooks/simulated_asteroid_ephem_rsp.ipynb +++ b/notebooks/sim_ephem/simulated_asteroid_ephem_rsp.ipynb @@ -20,11 +20,11 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:15.026373Z", - "iopub.status.busy": "2024-06-18T19:01:15.026079Z", - "iopub.status.idle": "2024-06-18T19:01:17.359920Z", - "shell.execute_reply": "2024-06-18T19:01:17.359236Z", - "shell.execute_reply.started": "2024-06-18T19:01:15.026354Z" + "iopub.execute_input": "2024-07-05T15:30:46.535448Z", + "iopub.status.busy": "2024-07-05T15:30:46.535291Z", + "iopub.status.idle": "2024-07-05T15:30:49.002871Z", + "shell.execute_reply": "2024-07-05T15:30:49.002322Z", + "shell.execute_reply.started": "2024-07-05T15:30:46.535434Z" } }, "outputs": [], @@ -49,11 +49,11 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:17.361390Z", - "iopub.status.busy": "2024-06-18T19:01:17.360891Z", - "iopub.status.idle": "2024-06-18T19:01:17.438804Z", - "shell.execute_reply": "2024-06-18T19:01:17.437930Z", - "shell.execute_reply.started": "2024-06-18T19:01:17.361368Z" + "iopub.execute_input": "2024-07-05T15:30:49.004023Z", + "iopub.status.busy": "2024-07-05T15:30:49.003641Z", + "iopub.status.idle": "2024-07-05T15:30:49.082021Z", + "shell.execute_reply": "2024-07-05T15:30:49.081508Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.004005Z" } }, "outputs": [], @@ -67,11 +67,11 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:17.439847Z", - "iopub.status.busy": "2024-06-18T19:01:17.439631Z", - "iopub.status.idle": "2024-06-18T19:01:17.442630Z", - "shell.execute_reply": "2024-06-18T19:01:17.442049Z", - "shell.execute_reply.started": "2024-06-18T19:01:17.439830Z" + "iopub.execute_input": "2024-07-05T15:30:49.082908Z", + "iopub.status.busy": "2024-07-05T15:30:49.082732Z", + "iopub.status.idle": "2024-07-05T15:30:49.085508Z", + "shell.execute_reply": "2024-07-05T15:30:49.085037Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.082894Z" } }, "outputs": [], @@ -85,11 +85,11 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:17.443664Z", - "iopub.status.busy": "2024-06-18T19:01:17.443428Z", - "iopub.status.idle": "2024-06-18T19:01:18.508689Z", - "shell.execute_reply": "2024-06-18T19:01:18.508031Z", - "shell.execute_reply.started": "2024-06-18T19:01:17.443646Z" + "iopub.execute_input": "2024-07-05T15:30:49.086275Z", + "iopub.status.busy": "2024-07-05T15:30:49.086098Z", + "iopub.status.idle": "2024-07-05T15:30:49.809760Z", + "shell.execute_reply": "2024-07-05T15:30:49.808044Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.086256Z" } }, "outputs": [], @@ -119,11 +119,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.509807Z", - "iopub.status.busy": "2024-06-18T19:01:18.509560Z", - "iopub.status.idle": "2024-06-18T19:01:18.661531Z", - "shell.execute_reply": "2024-06-18T19:01:18.660953Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.509790Z" + "iopub.status.busy": "2024-07-05T15:30:49.810323Z", + "iopub.status.idle": "2024-07-05T15:30:49.810571Z", + "shell.execute_reply": "2024-07-05T15:30:49.810457Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.810448Z" } }, "outputs": [], @@ -137,11 +136,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.663655Z", - "iopub.status.busy": "2024-06-18T19:01:18.663419Z", - "iopub.status.idle": "2024-06-18T19:01:18.683552Z", - "shell.execute_reply": "2024-06-18T19:01:18.682809Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.663638Z" + "iopub.status.busy": "2024-07-05T15:30:49.811051Z", + "iopub.status.idle": "2024-07-05T15:30:49.811253Z", + "shell.execute_reply": "2024-07-05T15:30:49.811163Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.811156Z" } }, "outputs": [], @@ -154,11 +152,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.685016Z", - "iopub.status.busy": "2024-06-18T19:01:18.684704Z", - "iopub.status.idle": "2024-06-18T19:01:18.699123Z", - "shell.execute_reply": "2024-06-18T19:01:18.698473Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.684988Z" + "iopub.status.busy": "2024-07-05T15:30:49.812059Z", + "iopub.status.idle": "2024-07-05T15:30:49.812280Z", + "shell.execute_reply": "2024-07-05T15:30:49.812190Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.812182Z" } }, "outputs": [], @@ -179,11 +176,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.700389Z", - "iopub.status.busy": "2024-06-18T19:01:18.700015Z", - "iopub.status.idle": "2024-06-18T19:01:18.715818Z", - "shell.execute_reply": "2024-06-18T19:01:18.715179Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.700367Z" + "iopub.status.busy": "2024-07-05T15:30:49.812815Z", + "iopub.status.idle": "2024-07-05T15:30:49.813013Z", + "shell.execute_reply": "2024-07-05T15:30:49.812926Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.812918Z" } }, "outputs": [], @@ -202,11 +198,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.716909Z", - "iopub.status.busy": "2024-06-18T19:01:18.716682Z", - "iopub.status.idle": "2024-06-18T19:01:18.743743Z", - "shell.execute_reply": "2024-06-18T19:01:18.743021Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.716892Z" + "iopub.status.busy": "2024-07-05T15:30:49.813642Z", + "iopub.status.idle": "2024-07-05T15:30:49.813847Z", + "shell.execute_reply": "2024-07-05T15:30:49.813758Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.813750Z" } }, "outputs": [], @@ -221,11 +216,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.745120Z", - "iopub.status.busy": "2024-06-18T19:01:18.744684Z", - "iopub.status.idle": "2024-06-18T19:01:18.761040Z", - "shell.execute_reply": "2024-06-18T19:01:18.760412Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.745100Z" + "iopub.status.busy": "2024-07-05T15:30:49.814629Z", + "iopub.status.idle": "2024-07-05T15:30:49.814844Z", + "shell.execute_reply": "2024-07-05T15:30:49.814751Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.814743Z" } }, "outputs": [], @@ -244,11 +238,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.762022Z", - "iopub.status.busy": "2024-06-18T19:01:18.761789Z", - "iopub.status.idle": "2024-06-18T19:01:18.775973Z", - "shell.execute_reply": "2024-06-18T19:01:18.775117Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.762000Z" + "iopub.status.busy": "2024-07-05T15:30:49.815488Z", + "iopub.status.idle": "2024-07-05T15:30:49.815705Z", + "shell.execute_reply": "2024-07-05T15:30:49.815616Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.815608Z" } }, "outputs": [], @@ -262,11 +255,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.777363Z", - "iopub.status.busy": "2024-06-18T19:01:18.776889Z", - "iopub.status.idle": "2024-06-18T19:01:18.796868Z", - "shell.execute_reply": "2024-06-18T19:01:18.796165Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.777340Z" + "iopub.status.busy": "2024-07-05T15:30:49.816210Z", + "iopub.status.idle": "2024-07-05T15:30:49.816409Z", + "shell.execute_reply": "2024-07-05T15:30:49.816322Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.816315Z" } }, "outputs": [], @@ -279,11 +271,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.798419Z", - "iopub.status.busy": "2024-06-18T19:01:18.797782Z", - "iopub.status.idle": "2024-06-18T19:01:18.816304Z", - "shell.execute_reply": "2024-06-18T19:01:18.815461Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.798397Z" + "iopub.status.busy": "2024-07-05T15:30:49.817161Z", + "iopub.status.idle": "2024-07-05T15:30:49.817360Z", + "shell.execute_reply": "2024-07-05T15:30:49.817271Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.817264Z" } }, "outputs": [], @@ -301,11 +292,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.817772Z", - "iopub.status.busy": "2024-06-18T19:01:18.817316Z", - "iopub.status.idle": "2024-06-18T19:01:18.833560Z", - "shell.execute_reply": "2024-06-18T19:01:18.832844Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.817749Z" + "iopub.status.busy": "2024-07-05T15:30:49.817935Z", + "iopub.status.idle": "2024-07-05T15:30:49.818130Z", + "shell.execute_reply": "2024-07-05T15:30:49.818043Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.818036Z" } }, "outputs": [], @@ -322,11 +312,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.834791Z", - "iopub.status.busy": "2024-06-18T19:01:18.834470Z", - "iopub.status.idle": "2024-06-18T19:01:18.848589Z", - "shell.execute_reply": "2024-06-18T19:01:18.847950Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.834772Z" + "iopub.status.busy": "2024-07-05T15:30:49.818889Z", + "iopub.status.idle": "2024-07-05T15:30:49.819093Z", + "shell.execute_reply": "2024-07-05T15:30:49.819005Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.818998Z" } }, "outputs": [], @@ -339,11 +328,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.849539Z", - "iopub.status.busy": "2024-06-18T19:01:18.849320Z", - "iopub.status.idle": "2024-06-18T19:01:18.865983Z", - "shell.execute_reply": "2024-06-18T19:01:18.865340Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.849513Z" + "iopub.status.busy": "2024-07-05T15:30:49.819774Z", + "iopub.status.idle": "2024-07-05T15:30:49.819972Z", + "shell.execute_reply": "2024-07-05T15:30:49.819883Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.819876Z" } }, "outputs": [], @@ -362,11 +350,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.867204Z", - "iopub.status.busy": "2024-06-18T19:01:18.866796Z", - "iopub.status.idle": "2024-06-18T19:01:18.880152Z", - "shell.execute_reply": "2024-06-18T19:01:18.879521Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.867185Z" + "iopub.status.busy": "2024-07-05T15:30:49.820507Z", + "iopub.status.idle": "2024-07-05T15:30:49.820733Z", + "shell.execute_reply": "2024-07-05T15:30:49.820637Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.820629Z" } }, "outputs": [], @@ -381,11 +368,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:18.881123Z", - "iopub.status.busy": "2024-06-18T19:01:18.880917Z", - "iopub.status.idle": "2024-06-18T19:01:51.982325Z", - "shell.execute_reply": "2024-06-18T19:01:51.981730Z", - "shell.execute_reply.started": "2024-06-18T19:01:18.881108Z" + "iopub.status.busy": "2024-07-05T15:30:49.821352Z", + "iopub.status.idle": "2024-07-05T15:30:49.821568Z", + "shell.execute_reply": "2024-07-05T15:30:49.821465Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.821457Z" } }, "outputs": [], @@ -418,11 +404,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:51.983360Z", - "iopub.status.busy": "2024-06-18T19:01:51.983151Z", - "iopub.status.idle": "2024-06-18T19:01:52.000157Z", - "shell.execute_reply": "2024-06-18T19:01:51.999473Z", - "shell.execute_reply.started": "2024-06-18T19:01:51.983345Z" + "iopub.status.busy": "2024-07-05T15:30:49.822161Z", + "iopub.status.idle": "2024-07-05T15:30:49.822363Z", + "shell.execute_reply": "2024-07-05T15:30:49.822271Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.822264Z" } }, "outputs": [], @@ -442,11 +427,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:52.001741Z", - "iopub.status.busy": "2024-06-18T19:01:52.001163Z", - "iopub.status.idle": "2024-06-18T19:01:52.010357Z", - "shell.execute_reply": "2024-06-18T19:01:52.009551Z", - "shell.execute_reply.started": "2024-06-18T19:01:52.001719Z" + "iopub.status.busy": "2024-07-05T15:30:49.822933Z", + "iopub.status.idle": "2024-07-05T15:30:49.823129Z", + "shell.execute_reply": "2024-07-05T15:30:49.823041Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.823034Z" } }, "outputs": [], @@ -467,11 +451,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:52.013233Z", - "iopub.status.busy": "2024-06-18T19:01:52.012889Z", - "iopub.status.idle": "2024-06-18T19:01:52.025365Z", - "shell.execute_reply": "2024-06-18T19:01:52.024660Z", - "shell.execute_reply.started": "2024-06-18T19:01:52.013210Z" + "iopub.status.busy": "2024-07-05T15:30:49.823966Z", + "iopub.status.idle": "2024-07-05T15:30:49.824165Z", + "shell.execute_reply": "2024-07-05T15:30:49.824078Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.824071Z" } }, "outputs": [], @@ -486,11 +469,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:52.026475Z", - "iopub.status.busy": "2024-06-18T19:01:52.026241Z", - "iopub.status.idle": "2024-06-18T19:01:52.095939Z", - "shell.execute_reply": "2024-06-18T19:01:52.095204Z", - "shell.execute_reply.started": "2024-06-18T19:01:52.026457Z" + "iopub.status.busy": "2024-07-05T15:30:49.824676Z", + "iopub.status.idle": "2024-07-05T15:30:49.824877Z", + "shell.execute_reply": "2024-07-05T15:30:49.824783Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.824776Z" } }, "outputs": [], @@ -505,11 +487,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:52.096956Z", - "iopub.status.busy": "2024-06-18T19:01:52.096748Z", - "iopub.status.idle": "2024-06-18T19:01:52.100730Z", - "shell.execute_reply": "2024-06-18T19:01:52.099911Z", - "shell.execute_reply.started": "2024-06-18T19:01:52.096940Z" + "iopub.status.busy": "2024-07-05T15:30:49.825499Z", + "iopub.status.idle": "2024-07-05T15:30:49.825709Z", + "shell.execute_reply": "2024-07-05T15:30:49.825623Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.825615Z" } }, "outputs": [], @@ -524,11 +505,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:52.101667Z", - "iopub.status.busy": "2024-06-18T19:01:52.101450Z", - "iopub.status.idle": "2024-06-18T19:01:52.118900Z", - "shell.execute_reply": "2024-06-18T19:01:52.118257Z", - "shell.execute_reply.started": "2024-06-18T19:01:52.101652Z" + "iopub.status.busy": "2024-07-05T15:30:49.826270Z", + "iopub.status.idle": "2024-07-05T15:30:49.826473Z", + "shell.execute_reply": "2024-07-05T15:30:49.826380Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.826373Z" } }, "outputs": [], @@ -542,11 +522,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:52.120101Z", - "iopub.status.busy": "2024-06-18T19:01:52.119878Z", - "iopub.status.idle": "2024-06-18T19:01:52.364967Z", - "shell.execute_reply": "2024-06-18T19:01:52.364185Z", - "shell.execute_reply.started": "2024-06-18T19:01:52.120086Z" + "iopub.status.busy": "2024-07-05T15:30:49.826974Z", + "iopub.status.idle": "2024-07-05T15:30:49.827178Z", + "shell.execute_reply": "2024-07-05T15:30:49.827091Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.827084Z" } }, "outputs": [], @@ -592,11 +571,10 @@ "execution_count": null, "metadata": { "execution": { - "iopub.execute_input": "2024-06-18T19:01:52.366083Z", - "iopub.status.busy": "2024-06-18T19:01:52.365860Z", - "iopub.status.idle": "2024-06-18T19:01:52.558412Z", - "shell.execute_reply": "2024-06-18T19:01:52.557808Z", - "shell.execute_reply.started": "2024-06-18T19:01:52.366066Z" + "iopub.status.busy": "2024-07-05T15:30:49.827679Z", + "iopub.status.idle": "2024-07-05T15:30:49.827879Z", + "shell.execute_reply": "2024-07-05T15:30:49.827788Z", + "shell.execute_reply.started": "2024-07-05T15:30:49.827780Z" } }, "outputs": [], diff --git a/notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb b/notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb new file mode 100644 index 0000000..72e6f1b --- /dev/null +++ b/notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb @@ -0,0 +1,910 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Testing DP0.3 ephemerides and coordinates\n", + "By Jamie Robinson" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook retrieves the observations and orbital parameters of an object in DP0.3. We compare the ephemerides and coordinates from DP0.3 to the output from propagating the orbit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:33:34.945521Z", + "iopub.status.busy": "2024-07-05T15:33:34.945104Z", + "iopub.status.idle": "2024-07-05T15:33:39.271163Z", + "shell.execute_reply": "2024-07-05T15:33:39.270462Z", + "shell.execute_reply.started": "2024-07-05T15:33:34.945496Z" + } + }, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "from astropy import units as u\n", + "from astropy.coordinates import SkyCoord, GCRS\n", + "from sbpy.data import Orbit\n", + "from astropy.time import Time\n", + "from astropy.coordinates import solar_system_ephemeris\n", + "from sbpy.photometry import HG\n", + "from astropy.table import QTable\n", + "\n", + "from lsst.rsp import get_tap_service\n", + "from astroquery.jplhorizons import Horizons" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:33:39.272663Z", + "iopub.status.busy": "2024-07-05T15:33:39.272094Z", + "iopub.status.idle": "2024-07-05T15:33:39.333026Z", + "shell.execute_reply": "2024-07-05T15:33:39.332159Z", + "shell.execute_reply.started": "2024-07-05T15:33:39.272639Z" + } + }, + "outputs": [], + "source": [ + "service = get_tap_service(\"ssotap\")\n", + "assert service is not None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:33:39.334212Z", + "iopub.status.busy": "2024-07-05T15:33:39.333990Z", + "iopub.status.idle": "2024-07-05T15:33:39.337089Z", + "shell.execute_reply": "2024-07-05T15:33:39.336551Z", + "shell.execute_reply.started": "2024-07-05T15:33:39.334195Z" + } + }, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "# ssoid = \"6098332225018\" # good test object\n", + "# ssoid = \"8268570668335894776\" # NEO\n", + "ssoid = \"-3649348068486548794\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-05T15:33:39.338054Z", + "iopub.status.busy": "2024-07-05T15:33:39.337859Z", + "iopub.status.idle": "2024-07-05T15:33:40.267531Z", + "shell.execute_reply": "2024-07-05T15:33:40.264933Z", + "shell.execute_reply.started": "2024-07-05T15:33:39.338038Z" + } + }, + "outputs": [], + "source": [ + "query = \"\"\"\n", + "SELECT\n", + " *\n", + "FROM\n", + " dp03_catalogs_10yr.DiaSource as dia\n", + "INNER JOIN\n", + " dp03_catalogs_10yr.SSSource as sss\n", + "ON\n", + " dia.diaSourceId = sss.diaSourceId\n", + "WHERE\n", + " dia.ssObjectId={}\n", + "ORDER by dia.ssObjectId\n", + "\"\"\".format(\n", + " ssoid\n", + ")\n", + "\n", + "df_obs = service.search(query).to_table().to_pandas()\n", + "df_obs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.268400Z", + "iopub.status.idle": "2024-07-05T15:33:40.268689Z", + "shell.execute_reply": "2024-07-05T15:33:40.268574Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.268562Z" + } + }, + "outputs": [], + "source": [ + "# calculate elongation angle\n", + "R = df_obs[\"heliocentricDist\"]\n", + "Delta = df_obs[\"topocentricDist\"]\n", + "alpha = np.radians(df_obs[\"phaseAngle\"])\n", + "\n", + "R_E = np.sqrt((R * R) + (Delta * Delta) - (2.0 * R * Delta * np.cos(alpha)))\n", + "df_obs[\"elong\"] = np.degrees(np.arccos(((R_E * R_E) + (Delta * Delta) - (R * R)) / (2.0 * R_E * Delta)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.269564Z", + "iopub.status.idle": "2024-07-05T15:33:40.269819Z", + "shell.execute_reply": "2024-07-05T15:33:40.269697Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.269688Z" + } + }, + "outputs": [], + "source": [ + "# put obs in time order\n", + "df_obs = df_obs.sort_values(\"midPointMjdTai\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.270771Z", + "iopub.status.idle": "2024-07-05T15:33:40.271049Z", + "shell.execute_reply": "2024-07-05T15:33:40.270945Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.270933Z" + } + }, + "outputs": [], + "source": [ + "# get orbit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.271857Z", + "iopub.status.idle": "2024-07-05T15:33:40.272093Z", + "shell.execute_reply": "2024-07-05T15:33:40.271983Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.271974Z" + } + }, + "outputs": [], + "source": [ + "results = service.search(\"SELECT * FROM dp03_catalogs_10yr.MPCORB \" \"WHERE ssObjectId={}\".format(ssoid))\n", + "df_orb = results.to_table().to_pandas()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.273422Z", + "iopub.status.idle": "2024-07-05T15:33:40.273664Z", + "shell.execute_reply": "2024-07-05T15:33:40.273559Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.273548Z" + } + }, + "outputs": [], + "source": [ + "df_orb" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Resample the sparse observations by propagating orbit with sbpy & oorb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.274352Z", + "iopub.status.idle": "2024-07-05T15:33:40.274586Z", + "shell.execute_reply": "2024-07-05T15:33:40.274472Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.274463Z" + } + }, + "outputs": [], + "source": [ + "# check orbit\n", + "e = df_orb[\"e\"]\n", + "incl = df_orb[\"incl\"]\n", + "q = df_orb[\"q\"]\n", + "a = q / (1.0 - e)\n", + "Q = a * (1.0 + e)\n", + "print(a, e, incl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.275180Z", + "iopub.status.idle": "2024-07-05T15:33:40.275404Z", + "shell.execute_reply": "2024-07-05T15:33:40.275301Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.275292Z" + } + }, + "outputs": [], + "source": [ + "df_orb[\"a\"] = a\n", + "df_orb[\"Q\"] = Q\n", + "df_orb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.276437Z", + "iopub.status.idle": "2024-07-05T15:33:40.276674Z", + "shell.execute_reply": "2024-07-05T15:33:40.276573Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.276563Z" + } + }, + "outputs": [], + "source": [ + "# Calculate some extra orbital elements\n", + "df_orb[\"P\"] = df_orb[\"a\"] ** (3.0 / 2.0) # orbital period in years\n", + "df_orb[\"n\"] = 360.0 / (df_orb[\"P\"] * 365.25) # mean motion in deg/day\n", + "df_orb[\"M\"] = (\n", + " df_orb[\"n\"] * (df_orb[\"epoch\"] - df_orb[\"tperi\"])\n", + ") % 360 # angles must be in correct range otherwise sbpy/pyoorb freak out\n", + "df_orb[[\"a\", \"P\", \"n\", \"M\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.277508Z", + "iopub.status.idle": "2024-07-05T15:33:40.277727Z", + "shell.execute_reply": "2024-07-05T15:33:40.277634Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.277626Z" + } + }, + "outputs": [], + "source": [ + "# rename columns for consistency\n", + "df_orb = df_orb.rename(columns={\"node\": \"Omega\", \"peri\": \"w\"})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.278432Z", + "iopub.status.idle": "2024-07-05T15:33:40.278645Z", + "shell.execute_reply": "2024-07-05T15:33:40.278549Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.278541Z" + } + }, + "outputs": [], + "source": [ + "df_orb[[\"a\", \"e\", \"incl\", \"Omega\", \"w\", \"M\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.289637Z", + "iopub.status.idle": "2024-07-05T15:33:40.290438Z", + "shell.execute_reply": "2024-07-05T15:33:40.290280Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.290263Z" + } + }, + "outputs": [], + "source": [ + "# # create an sbpy oorb object from dataframe via QTable\n", + "# tab = QTable.from_pandas(\n", + "# df_orb[[\"a\", \"e\", \"incl\", \"Omega\", \"w\", \"M\"]],\n", + "# units={\"a\": u.au, \"incl\": u.deg, \"Omega\": u.deg, \"w\": u.deg, \"M\": u.deg},\n", + "# )\n", + "# orbit = Orbit.from_table(tab)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.291224Z", + "iopub.status.idle": "2024-07-05T15:33:40.291730Z", + "shell.execute_reply": "2024-07-05T15:33:40.291609Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.291597Z" + } + }, + "outputs": [], + "source": [ + "# # oorb requires certain extra fields\n", + "# orbit[\"epoch\"] = Time(Time(df_orb[\"epoch\"], format=\"mjd\").jd, format=\"jd\")\n", + "# orbit[\"targetname\"] = np.array(df_orb[\"ssObjectId\"]).astype(str)\n", + "# orbit[\"H\"] = df_orb[\"mpcH\"] * u.mag\n", + "# orbit[\"G\"] = df_orb[\"mpcG\"] * u.dimensionless_unscaled" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.292524Z", + "iopub.status.idle": "2024-07-05T15:33:40.292763Z", + "shell.execute_reply": "2024-07-05T15:33:40.292661Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.292652Z" + } + }, + "outputs": [], + "source": [ + "# orbit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.293660Z", + "iopub.status.idle": "2024-07-05T15:33:40.293932Z", + "shell.execute_reply": "2024-07-05T15:33:40.293825Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.293812Z" + } + }, + "outputs": [], + "source": [ + "# # define a set of JD times to propagate the orbital elements to\n", + "# N = 1000\n", + "# times = Time(\n", + "# Time(np.linspace(np.amin(df[\"midPointMjdTai\"]), np.amax(df[\"midPointMjdTai\"]), N), format=\"mjd\").jd,\n", + "# format=\"jd\",\n", + "# )\n", + "# times[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.294706Z", + "iopub.status.idle": "2024-07-05T15:33:40.294965Z", + "shell.execute_reply": "2024-07-05T15:33:40.294855Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.294845Z" + } + }, + "outputs": [], + "source": [ + "# # create an empty dataframe to hold resampled observations\n", + "# df_dense = pd.DataFrame()\n", + "# df_dense[\"midPointMjdTai\"] = times.mjd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.296032Z", + "iopub.status.idle": "2024-07-05T15:33:40.296273Z", + "shell.execute_reply": "2024-07-05T15:33:40.296163Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.296154Z" + } + }, + "outputs": [], + "source": [ + "# # propagate the orbit forward in time.\n", + "# # probably a better way to do this but I can't get oo_propagate to work with a Time list right now\n", + "# # see: https://github.com/NASA-Planetary-Science/sbpy/issues/341\n", + "\n", + "# df_pos = pd.DataFrame() # empty dataframe to hold cartesian coordinates\n", + "\n", + "# for i in range(len(times)):\n", + "# print(i)\n", + "# prop_elem = orbit.oo_propagate(times[i]) # propagate the orbit to the selected time step\n", + "# del prop_elem.table[\n", + "# \"orbtype\"\n", + "# ] # orbtype is added as int, sbpy freaks out so delete the orbtype and then _to_oo works it out\n", + "# print(\"propagate\")\n", + "# statevec = prop_elem.oo_transform(\"CART\") # transform from orbital elements to cartesian\n", + "# print(\"transform\")\n", + "\n", + "# # append new cartesian coordinates to the dataframe\n", + "# _df_statevec = statevec.table.to_pandas()\n", + "# df_pos = pd.concat((df_pos, _df_statevec))\n", + "\n", + "# df_pos.reset_index(drop=True, inplace=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.296898Z", + "iopub.status.idle": "2024-07-05T15:33:40.297138Z", + "shell.execute_reply": "2024-07-05T15:33:40.297029Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.297018Z" + } + }, + "outputs": [], + "source": [ + "# df_pos" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Use astropy coordinates to transform between all the coordinate systems" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.298026Z", + "iopub.status.idle": "2024-07-05T15:33:40.298250Z", + "shell.execute_reply": "2024-07-05T15:33:40.298152Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.298143Z" + } + }, + "outputs": [], + "source": [ + "# # define heliocentric cartesian coordinates\n", + "# c_xyz_hel = SkyCoord(\n", + "# x=np.array(df_pos[\"x\"]),\n", + "# y=np.array(df_pos[\"y\"]),\n", + "# z=np.array(df_pos[\"z\"]),\n", + "# unit=\"AU\",\n", + "# representation_type=\"cartesian\",\n", + "# frame=\"heliocentrictrueecliptic\",\n", + "# )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.298942Z", + "iopub.status.idle": "2024-07-05T15:33:40.299172Z", + "shell.execute_reply": "2024-07-05T15:33:40.299073Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.299064Z" + } + }, + "outputs": [], + "source": [ + "# # transform to heliocentric ecliptic coords\n", + "# c_ecl_hel = c_xyz_hel.copy()\n", + "# c_ecl_hel.representation_type = \"spherical\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.299978Z", + "iopub.status.idle": "2024-07-05T15:33:40.300217Z", + "shell.execute_reply": "2024-07-05T15:33:40.300117Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.300108Z" + } + }, + "outputs": [], + "source": [ + "# # transform to geocentric equatorial coords (times required to calculate Earth position)\n", + "# with solar_system_ephemeris.set(\"jpl\"):\n", + "# c_eq_geo = c_xyz_hel.transform_to(GCRS(obstime=times))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.301111Z", + "iopub.status.idle": "2024-07-05T15:33:40.301359Z", + "shell.execute_reply": "2024-07-05T15:33:40.301243Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.301233Z" + } + }, + "outputs": [], + "source": [ + "# # transform to geocentric cartesian coords\n", + "# c_xyz_geo = c_eq_geo.copy()\n", + "# c_xyz_geo.representation_type = \"cartesian\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.301927Z", + "iopub.status.idle": "2024-07-05T15:33:40.302147Z", + "shell.execute_reply": "2024-07-05T15:33:40.302042Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.302034Z" + } + }, + "outputs": [], + "source": [ + "# # transform from geo equatorial (ra, dec) to geo ecliptic (lon, lat)\n", + "# c_ecl_geo = c_eq_geo.transform_to(\"geocentrictrueecliptic\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.302981Z", + "iopub.status.idle": "2024-07-05T15:33:40.303210Z", + "shell.execute_reply": "2024-07-05T15:33:40.303111Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.303102Z" + } + }, + "outputs": [], + "source": [ + "# # plot the propagated cartesian positions against the database values\n", + "\n", + "# x_plot = \"midPointMjdTai\"\n", + "# y_plot1 = \"heliocentricX\"\n", + "# y_plot2 = \"heliocentricY\"\n", + "# y_plot3 = \"heliocentricZ\"\n", + "# df_plot = df\n", + "\n", + "# fig = plt.figure()\n", + "# gs = gridspec.GridSpec(1, 1)\n", + "# ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "# x = ax1.scatter(df_plot[x_plot], df_plot[y_plot1], label=y_plot1)\n", + "# x = ax1.scatter(df_plot[x_plot], df_plot[y_plot2], label=y_plot2)\n", + "# x = ax1.scatter(df_plot[x_plot], df_plot[y_plot3], label=y_plot3)\n", + "\n", + "# ax1.plot(times.mjd, c_xyz_hel.x)\n", + "# ax1.plot(times.mjd, c_xyz_hel.y)\n", + "# ax1.plot(times.mjd, c_xyz_hel.z)\n", + "\n", + "# ax1.set_xlabel(x_plot)\n", + "# ax1.set_ylabel(\"distance\")\n", + "# ax1.legend()\n", + "\n", + "# plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are some deviations of x and y positions, probably due to slightly different reference frames and methods of propagating orbits.\n", + "\n", + "There is something wrong with database z positions, will be fixed soon!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.303807Z", + "iopub.status.idle": "2024-07-05T15:33:40.304051Z", + "shell.execute_reply": "2024-07-05T15:33:40.303949Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.303940Z" + } + }, + "outputs": [], + "source": [ + "# # the ecliptic coordinates look good!\n", + "\n", + "# x_plot = \"midPointMjdTai\"\n", + "# y_plot1 = \"eclipticLambda\"\n", + "# y_plot2 = \"eclipticBeta\"\n", + "# df_plot = df\n", + "\n", + "# fig = plt.figure()\n", + "# gs = gridspec.GridSpec(1, 1)\n", + "# ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "# x = ax1.scatter(df_plot[x_plot], df_plot[y_plot1], label=y_plot1)\n", + "# x = ax1.scatter(df_plot[x_plot], df_plot[y_plot2], label=y_plot2)\n", + "\n", + "# ax1.plot(times.mjd, c_ecl_geo.lon.degree)\n", + "# ax1.plot(times.mjd, c_ecl_geo.lat.degree)\n", + "\n", + "# ax1.set_xlabel(x_plot)\n", + "# ax1.set_ylabel(\"angle\")\n", + "# ax1.legend()\n", + "\n", + "# plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.304969Z", + "iopub.status.idle": "2024-07-05T15:33:40.305201Z", + "shell.execute_reply": "2024-07-05T15:33:40.305099Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.305090Z" + } + }, + "outputs": [], + "source": [ + "# query JPL" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.305853Z", + "iopub.status.idle": "2024-07-05T15:33:40.306082Z", + "shell.execute_reply": "2024-07-05T15:33:40.305973Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.305966Z" + } + }, + "outputs": [], + "source": [ + "target = df_orb.iloc[0][\"fullDesignation\"].split(\"2011 \")[-1]\n", + "print(target)\n", + "\n", + "site = \"X05\" # Roques de los Muchachos\n", + "times = {\"start\": \"2023-10-01 04:00\", \"stop\": \"2033-10-01 04:00\", \"step\": \"1day\"} # dates to query\n", + "obj = Horizons(id=target, location=site, epochs=times)\n", + "eph = obj.ephemerides()\n", + "df_eph = eph.to_pandas()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.307034Z", + "iopub.status.idle": "2024-07-05T15:33:40.307249Z", + "shell.execute_reply": "2024-07-05T15:33:40.307156Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.307148Z" + } + }, + "outputs": [], + "source": [ + "obj = Horizons(id=target, epochs=times, location=\"399\")\n", + "vec = obj.vectors()\n", + "df_vec_earth = vec.to_pandas()\n", + "\n", + "obj = Horizons(id=target, epochs=times, location=\"@10\")\n", + "vec = obj.vectors()\n", + "df_vec_sun = vec.to_pandas()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.308110Z", + "iopub.status.idle": "2024-07-05T15:33:40.308328Z", + "shell.execute_reply": "2024-07-05T15:33:40.308227Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.308220Z" + } + }, + "outputs": [], + "source": [ + "df_vec_earth = df_vec_earth.rename({\"x\": \"topocentricX\", \"y\": \"topocentricY\", \"z\": \"topocentricZ\"}, axis=1)\n", + "df_vec_sun = df_vec_sun.rename({\"x\": \"heliocentricX\", \"y\": \"heliocentricY\", \"z\": \"heliocentricZ\"}, axis=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.309065Z", + "iopub.status.idle": "2024-07-05T15:33:40.309290Z", + "shell.execute_reply": "2024-07-05T15:33:40.309190Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.309182Z" + } + }, + "outputs": [], + "source": [ + "df_vec = df_vec_earth[\n", + " [\"targetname\", \"datetime_jd\", \"datetime_str\", \"topocentricX\", \"topocentricY\", \"topocentricZ\"]\n", + "].merge(\n", + " df_vec_sun[\n", + " [\"targetname\", \"datetime_jd\", \"datetime_str\", \"heliocentricX\", \"heliocentricY\", \"heliocentricZ\"]\n", + " ],\n", + " on=[\"targetname\", \"datetime_jd\", \"datetime_str\"],\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.310032Z", + "iopub.status.idle": "2024-07-05T15:33:40.310264Z", + "shell.execute_reply": "2024-07-05T15:33:40.310165Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.310156Z" + } + }, + "outputs": [], + "source": [ + "df_eph[\"datetime_mjd\"] = Time(df_eph[\"datetime_jd\"], format=\"jd\").mjd\n", + "df_vec[\"datetime_mjd\"] = Time(df_vec[\"datetime_jd\"], format=\"jd\").mjd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.311155Z", + "iopub.status.idle": "2024-07-05T15:33:40.311405Z", + "shell.execute_reply": "2024-07-05T15:33:40.311286Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.311277Z" + } + }, + "outputs": [], + "source": [ + "night_mask = (~np.isin(df_eph[\"solar_presence\"], [\"*\", \"N\", \"C\"])) & (df_eph[\"EL\"] > 0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.312011Z", + "iopub.status.idle": "2024-07-05T15:33:40.312264Z", + "shell.execute_reply": "2024-07-05T15:33:40.312149Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.312141Z" + } + }, + "outputs": [], + "source": [ + "x_plot = \"midPointMjdTai\"\n", + "y_plot = \"elong\"\n", + "df_plot = df_obs\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(df_plot[x_plot], df_plot[y_plot], label=\"DP0.3\")\n", + "ax1.plot(df_eph[\"datetime_mjd\"], df_eph[\"elong\"], c=\"r\", label=\"JPL\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.313021Z", + "iopub.status.idle": "2024-07-05T15:33:40.313245Z", + "shell.execute_reply": "2024-07-05T15:33:40.313148Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.313139Z" + } + }, + "outputs": [], + "source": [ + "x_plot = \"midPointMjdTai\"\n", + "y_plot = \"phaseAngle\"\n", + "df_plot = df_obs\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.scatter(df_plot[x_plot], df_plot[y_plot], label=\"DP0.3\")\n", + "ax1.plot(df_eph[\"datetime_mjd\"], df_eph[\"alpha\"], c=\"r\", label=\"JPL\")\n", + "ax1.scatter(df_eph[night_mask][\"datetime_mjd\"], df_eph[night_mask][\"alpha\"], edgecolor=\"C1\", facecolor=\"none\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.status.busy": "2024-07-05T15:33:40.314008Z", + "iopub.status.idle": "2024-07-05T15:33:40.314240Z", + "shell.execute_reply": "2024-07-05T15:33:40.314140Z", + "shell.execute_reply.started": "2024-07-05T15:33:40.314130Z" + } + }, + "outputs": [], + "source": [ + "x_plot = \"heliocentricX\"\n", + "y_plot = \"heliocentricY\"\n", + "c_plot = \"midPointMjdTai\"\n", + "df_plot = df_obs\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "s1 = ax1.scatter(df_plot[x_plot], df_plot[y_plot], c=df_plot[c_plot])\n", + "cbar1 = plt.colorbar(s1)\n", + "\n", + "ax1.scatter(0, 0, marker=\"+\", c=\"k\")\n", + "circle1 = plt.Circle((0, 0), 1.0, edgecolor=\"r\", facecolor=\"none\", label=\"1au\")\n", + "ax1.add_patch(circle1)\n", + "\n", + "mask = df_plot[\"phaseAngle\"] > 120\n", + "_df_plot = df_plot[mask]\n", + "ax1.scatter(_df_plot[x_plot], _df_plot[y_plot], facecolor=\"none\", edgecolor=\"r\")\n", + "\n", + "ax1.set_aspect(\"equal\")\n", + "ax1.legend()\n", + "\n", + "cbar1.set_label(c_plot)\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "lsst-sbpy", + "language": "python", + "name": "lsst-sbpy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 5a8d15cc058962a943af927d6b64a5e5b0e04505 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Mon, 8 Jul 2024 13:56:54 +0000 Subject: [PATCH 8/8] update blurb --- notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb b/notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb index 72e6f1b..fb45ddb 100644 --- a/notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb +++ b/notebooks/sim_ephem/simulated_asteroid_ephem_rsp_horizons.ipynb @@ -12,7 +12,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This notebook retrieves the observations and orbital parameters of an object in DP0.3. We compare the ephemerides and coordinates from DP0.3 to the output from propagating the orbit." + "This notebook retrieves the observations and orbital parameters of an object in DP0.3. We compare the ephemerides and coordinates from DP0.3 to the JPL Horizons values." ] }, {