Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Palomar Plate Tutorial #108

Merged
merged 2 commits into from
Sep 1, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

### Added

- Added `kete.RectangleFOV.from_wcs`, allowing the construction of a FOV from a given
Astropy WCS object.
- Added `kete.conversion.bin_data`, which allows for binning matrix data such as images.
- Added tutorial showing precovery of an asteroid from a 1950's glass plate observation
done at the Palomar Observatory.
- Added sunshield rotation calculation for NEO Surveyor.


Expand Down
Binary file added docs/data/full_frame.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/data/full_frame_annotated.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/data/full_frame_annotated_zoom.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
1 change: 1 addition & 0 deletions docs/tutorials/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -10,5 +10,6 @@ These are more in-depth worked examples.
background
propagation
neatm
palomar
kona
wise
224 changes: 224 additions & 0 deletions docs/tutorials/palomar.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,224 @@
Palomar 1950 Plate
==================

Identification of an asteroid in a glass plate exposure from 1950 Palomar Mountain.

Background
----------

On December 6th and 8th 1950, Rudolph Minkowski captured 2 long exposure images
onto glass plates of the galaxy NGC2403 looking for time varying sources. On the
cover slip for these plates was written "Asteroid", and an arrow was drawn on the
plate identifying the streaking asteroid on the plates themselves. Details of the
observations are included in
`Tammann & Sandage (1968) <https://doi.org/10.1086/149487>`_ in their Table A1,
although plate PH-98M is incorrectly listed as PH-98H (with the ending letter
indicating the observer). These plates were exposed for approximately 30 minutes
(although cloud cover reduced the exposure toward the end for one plate).

These plates are stored at the Carnegie Science Observatories in Pasadena in their
archives. One of Carnegie's Archivists, Kit Whitten, provided a scanned copy of
the frames to Joe Masiero at CalTech IPAC around 2018. Masiero provided them for
the following analysis.

Here we will identify that asteroid which Minkowski captured over 70 years ago.


Pre-processing
--------------

Some pre-processing was done on the scanned image before this analysis. First the
data was binned by a factor of 8 in both axis to reduce the file size. Second it
was passed to astrometry.net for a plate solution. This gives a valid WCS to work
with, which is practically a requirement. Since the image was as scan of a
film negative, it was also inverted.

Limitations
-----------

- Astrometry.net is using the current star catalog, meaning if there are stars in
the field which have moved measurably, then the astrometry will be imperfect.
- Observer position before 1972 is somewhat tricky to compute, and is currently not
well supported in Kete. This is a result of NASA NAIF (Navigation and Ancillary
Information Facility) only producing Earth orientation files back to about 1972.
This information can be calculated by other means with less precision, but this
is not currently implemented. As a result of this, the center of the Earth is used
as an approximation for the observers position, but this results in about an
arc-minute offset in calculated positions. Adding a better approximation of the
observer location shifts the results to sub PSF accuracy.
- Record keeping of the exact time appears to be off by a few minutes, though this
may be an artifact of the observer position being approximate.


FITs File
---------

First we load a FITs file, and grab frame information from its header.

.. code-block:: python

import kete
import matplotlib.pyplot as plt
import numpy as np

from astropy.io import fits
from astropy.wcs import WCS


# load the fits file and associated WCS
frame = fits.open("data/19501206-NGC2403-reduced.fits.gz")[0]
wcs = WCS(frame.header)

# Time of observation, approximately 3 am local time
# Taken from table A1 of the paper cited above.
jd = kete.Time(2433622.977, scaling='utc').jd

earth = kete.spice.get_state("Earth", jd)
fov = kete.fov.RectangleFOV.from_wcs(wcs, earth)

# Plot the frame
kete.irsa.plot_fits_image(frame)
plt.title("Dec 6 1950 - Palomar Mountain")
plt.tight_layout()
plt.savefig("data/full_frame.png")

.. image:: ../data/full_frame.png
:alt: Raw frame scanned from the Dec 6 1950 plate.


MPC Orbit Data
--------------

Next we must collect orbit information from the Minor Planet Center (MPC).
We will load all known objects from their database, convert them to a kete State,
and propagate those states to the epoch near the FITs file epoch we opened.

We also trim out more distant and eccentric objects.

.. code-block:: python

table = kete.mpc.fetch_known_orbit_data()
table = table[[str(t).isdigit() for t in table['desig']]]
table = table[table['ecc'] < 0.5]
table = table[table['peri_dist'] < 4]
states = kete.mpc.table_to_states(table)


Quick Cut
---------

As we are only analyzing a single frame, and the speed of the object suggests it
is not close to Earth, we can apply some approximations to remove many objects from
contention.

.. note::

This step below can be completely skipped, it is only an optimization to make the
orbit propagation faster. This step only really applies for off ecliptic field of
views, as any FOV on the ecliptic will already have most of the main belt as
possible objects.

If we assume most objects stay in their orbital plane, but we simply don't know where
they are along the plane, we can use the two body approximation to compute the plane's
on sky position from the observers location.

Here we assume that the object is likely a main belt, or main belt adjacent object,
and in doing so we can assume an orbital period of between about 800 and 2500 days.

If we sample the orbit 360 times around the period of 2500 days, objects with that period
will be sampled at about 1 degree steps on the sky. We then check at each of these points
how close the object came to the center of the FOV.

.. code-block:: python

# Define a convenience function to see how far the objects are from the FOV
def cur_angles(states, fov):
"""
Given a list of states, compute how far they are from the center of a FOV
"""
pointing = fov.pointing
angles = []
obs_pos = fov.observer.pos
for state in states:
vec = state.pos - obs_pos
angles.append(pointing.angle_between(vec))
return np.array(angles)

n_steps = 360
period = 2500

jd = np.median([s.jd for s in states])
states = kete.propagate_two_body(states, jd)
best = cur_angles(states, fov)

for dt in np.linspace(0, 1, n_steps) * period:
states = kete.propagate_two_body(states, jd + dt)
cur = cur_angles(states, fov)
best = np.min([best, cur], axis=0)

mask = best < 3
print(f"Total of {sum(mask)} asteroids less than 3 degrees from the FOV")
states = kete.mpc.table_to_states(table[mask])

::

Total of 1148 asteroids less than 3 degrees from the FOV

Exact Calculation
-----------------

Compute the visible states, including large main belt asteroids in the computation.
For 1k objects this should take about a second.

.. code-block:: python

vis = kete.fov_state_check(states, [fov], include_asteroids=True)[0]

Results
-------

Plot the overlap of visible objects on the original frame.

Two objects are possible, 382632 and 151109, however the second is off the edge of
the glass plate itself so not visible.

.. code-block:: python

wcs = kete.irsa.plot_fits_image(frame, percentiles=(40, 99))
for idx in range(len(vis)):
vec = vis.obs_vecs()[idx]
kete.irsa.annotate_plot(wcs, vec, style='o', px_gap=10, text=vis[idx].desig)

plt.title("Dec 6 1950 - Annotated")
plt.savefig("data/full_frame_annotated.png")


.. image:: ../data/full_frame_annotated.png
:alt: Frame with both objects annotated.


Zoomed
------

Zooming in to the expected position of 382632, we see the original hand drawn arrow,
along with correctly labeled streaking asteroid. Note the slight offset, which from
analysis not performed here can be shown to be a result of our assumption that
the observer is at the center of the Earth. A better assumption about observer
position causes the alignment to match within the width of the blur.

.. code-block:: python

wcs = kete.irsa.plot_fits_image(frame, percentiles=(30, 99.9))
kete.irsa.annotate_plot(wcs,
vis.obs_vecs()[1],
style='L',
px_gap=5,
length=10,
text=" " + vis[1].desig)
kete.irsa.zoom_plot(wcs, vis.obs_vecs()[1])
plt.title(f"Dec 6 1950 - Annotated Zoom");
plt.savefig("data/full_frame_annotated_zoom.png")
plt.close()

.. image:: ../data/full_frame_annotated_zoom.png
:alt: Zoom in of the expected object.
1 change: 1 addition & 0 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ docs = ["sphinx",
"sphinx-gallery",
"sphinx_rtd_theme",
"click",
"pytz",
]

# Options for pytest
Expand Down
49 changes: 49 additions & 0 deletions src/kete/conversion.py
Original file line number Diff line number Diff line change
Expand Up @@ -195,6 +195,55 @@ def compute_tisserand(
)


def bin_data(data, bin_size=2, method="mean"):
"""
Bin data, reducing the size of a matrix.

Parameters
----------
data:
2 dimensional array of data to be reduced.
bin_size:
Binning size, binning of data is applied to both directions of the data.
If the bin size is not a multiple of the shape of the data, the data is
trimmed until it is.
method:
Which function is used to produce a single final value from the bins.
This includes 'mean', 'median', 'min', 'max'. NANs are ignored.
"""
funcs = {
"mean": np.nanmean,
"median": np.nanmedian,
"min": np.nanmin,
"max": np.nanmax,
}

# dont trust the user, check inputs
if bin_size <= 0 or int(bin_size) != bin_size:
raise ValueError("Bin size must be larger than 0 and an integer.")
bin_size = int(bin_size)
if method not in funcs:
raise ValueError("'method' must be one of: ", list(funcs.keys()))

data = np.asarray(data)
if data.ndim != 2:
raise ValueError("Array must be 2 dimensional.")

xdim, ydim = data.shape
if xdim % bin_size != 0:
xdim = (xdim // bin_size) * bin_size
logger.warning("x dimension is not a multiple of bin size, trimming to size.")
data = data[:xdim, :]
if ydim % bin_size != 0:
ydim = (ydim // bin_size) * bin_size
logger.warning("y dimension is not a multiple of bin size, trimming to size.")
data = data[:, :ydim]

reshaped = data.reshape(xdim // bin_size, bin_size, ydim // bin_size, bin_size)

return funcs[method](funcs[method](reshaped, axis=1), axis=2)


def dec_degrees_to_dms(dec: NDArray) -> Union[str, list[str]]:
"""
Convert a declination in degrees to a "degrees arcminutes arcseconds" string.
Expand Down
29 changes: 28 additions & 1 deletion src/kete/fov.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
fov_state_check,
fov_spk_check as _fov_spk_check,
)
from .vector import State
from .vector import State, Vector
from . import spice


Expand Down Expand Up @@ -55,3 +55,30 @@ def fov_spice_check(desigs: list[str], fovs) -> list[State]:
"""
obj_ids = [spice.name_lookup(n)[1] for n in desigs]
return _fov_spk_check(obj_ids, fovs)


def _from_wcs(wcs, obs: State) -> RectangleFOV:
"""
Construct a RectangleFOV from an astropy WCS along with an observer state.

parameters
----------
wcs:
An astropy WCS, this must include the shape of the array.
obs:
The observer position
"""
types = wcs.axis_type_names
if types == ["RA", "DEC"]:
x, y = wcs.array_shape
corners = wcs.pixel_to_world_values([0, 0, x, x], [0, y, y, 0])
vecs = [Vector.from_ra_dec(ra, dec) for ra, dec in zip(*corners)]
else:
raise NotImplementedError(
f"Support for WCS with frame {types} is not currently supported."
)
return RectangleFOV.from_corners(vecs, obs)


# Monkey patch the rust class with a new constructor.
RectangleFOV.from_wcs = staticmethod(_from_wcs)
7 changes: 4 additions & 3 deletions src/kete/rust/propagation.rs
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
use itertools::Itertools;
use kete_core::{
errors::Error,
propagation::{self, moid, NonGravModel},
spice::{self, get_spk_singleton},
state::State,
time::{scales::TDB, Time},
};
use itertools::Itertools;
use pyo3::{pyfunction, PyResult, Python};
use rayon::prelude::*;

Expand Down Expand Up @@ -165,8 +165,9 @@ pub fn propagation_n_body_spk_py(
///
/// This returns two lists of states:
/// - First one contains the states of the objects at the end of the integration
/// - Second contains the states of the planets at the end of the integration, which may
/// be used as input for continuing the integration.
/// - Second contains the states of the planets at the end of the integration.
///
/// The second set of states may be used as input for continuing the integration.
///
/// Parameters
/// ----------
Expand Down