ortk95 / planetmapper

PlanetMapper: An open source Python package for visualising, navigating and mapping Solar System observations
https://planetmapper.readthedocs.io
MIT License
10 stars 1 forks source link

Issue with plotting when target at high declinations #323

Closed goatsneez closed 7 months ago

goatsneez commented 7 months ago

This is a wonderful package. I had to work on a similar topic, but this could be a time saver if it can be used to plot also targets close to the observer.

There seems to be an issue of plotting wire_frame figures when the observer is closer to the target. In this example it is the JUICE spacecraft observing Jupiter. (It works fine for times when the target is rather far-away maybe 50x planet diameter or so -- but that just educate guess).

Minimal code for reproducing this bug is pasted below, it relies on the JUICE spice-kernels:

mkpath = '../../../kernels/juice/kernels/mk/juice_plan.tm' plm.base.prevent_kernel_loading() spice.furnsh(mkpath) body_jup = plm.Body('Jupiter','2032-09-24T05:30:00',\ observer='JUICE_SPACECRAFT',\ observer_frame='JUICE_SPACECRAFT_PLAN')

ortk95 commented 7 months ago

Thanks for the bug report! I can reproduce the issue, so there's definitely something going wrong - I'll see if I can work out what the issue is...

ortk95 commented 7 months ago

Code to reproduce issue

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
import spiceypy as spice

import planetmapper

mkpath = './scratchpad/JUICE/kernels/mk/juice_plan.tm'
planetmapper.base.prevent_kernel_loading()
spice.kclear()
spice.furnsh(mkpath)

times = [
    '2032-09-01T00:00:00',
    '2032-09-24T05:30:00',
    '2032-09-24T06:30:00',
]

for time in times:
    body = planetmapper.Body(
        'Jupiter',
        time,
        observer='JUICE_SPACECRAFT',
        observer_frame='JUICE_SPACECRAFT_PLAN',
    )
    body.plot_wireframe_km()
    plt.show()

image image image

ortk95 commented 7 months ago

Looking at the RA/Dec plots, it looks like the issue is caused by Jupiter being at a declination close 90, where the coordinate singularity causes significant warping of the RA/Dec coordinates. This is probably always going to happen somewhere when plotting spherical coordinates on a cartesian axis, but we haven't encountered it with real data before, as most long distance telescope observations are close to the ecliptic so will have low declination angles.

image

The singularity is inherent in the radec coordinates so may not be be fixable (without breaking the simplicity of the plots for other datasets at least). However, it should be possible to change how the km/xy coordinate systems are calculated, to ensure they always work, even at high declination angles.

Things to do

ortk95 commented 7 months ago

@goatsneez for your specific example, the issue only occurs when using the JUICE_SPACECRAFT_PLAN frame. For the (default) J2000 frame, the declinations are much smaller which avoids the singularity at the celestial pole. So unless you specifically need to be using the JUICE_SPACECRAFT_PLAN frame, you may be okay just using J2000 for now as a temporary workaround while this problem is fixed:

image

goatsneez commented 7 months ago

Thanks for the quick analysis; glad that the issue is understood, and good that there is an some work around.

As for my own needs, I would need the S/C frame since I would like to plot other pointings on that plot (instrumental FOVs, etc). The J2000 frame would be more cumbersome to look at in that context.

ortk95 commented 7 months ago

Hi @goatsneez! I've just released PlanetMapper 1.9.0 which includes some updates that may be able to help with your specific issue.

There's a new relative angular coordinate system, which is very similar to the RA/Dec coordinates, but with a different origin for the spherical coordinate system. That means that you can avoid the coordinate singularity at the celestial poles by choosing an origin where the target isn't near the pole. The relevant bits of the documentation are probably Body.radec2angular and Body.plot_wireframe_angular, and some small updates to the examples page. I've also added this problem to the common issues page, so thanks very much for bringing it to my attention!

By default, the angular coordinates are centred on the target body, so it should never appear distorted, but you could also specify some custom origin if you want the coordinates to be relative to the spacecraft frame. If I remember right, MAJIS, JANUS, UVS etc. all have a similar orientation, so there's probably some spacecraft axis that could make sense as an origin for the relative angular coordinates.

You can then use the Body.radec2angular function to convert RA/Dec coordinates (defining the instrument FOV) to the angular coordinates that you can plot with the wireframe. For example, here's a very quick mockup of how it would work using the different coordinate systems:

import matplotlib.pyplot as plt

import planetmapper

body = planetmapper.Body('jupiter', '2024-01-01')

# Define an example "instrument" FOV as a rectangle in RA/Dec coordinates
ra0, dec0 = 33.33, 12.15  # RA/Dec corners of the FOV
ra1, dec1 = 33.35, 12.16
ras = [ra0, ra1, ra1, ra0, ra0]  # Arrays to plot the FOV rectangle
decs = [dec0, dec0, dec1, dec1, dec0]

# Plot in RA/Dec frame
# - Simple, but the target may appear distorted at very high declinations.
ax = body.plot_wireframe_radec()
ax.plot(ras, decs, color='r')
plt.show()

# Plot in angular frame centred on target
# - The target should always appear undistorted.
ax = body.plot_wireframe_angular()
x, y = zip(*[body.radec2angular(ra, dec) for ra, dec in zip(ras, decs)])
ax.plot(x, y, color='r')
plt.show()

# Plot in angular frame centred on some specific RA/Dec (e.g. an instrument pointing)
# - This will be undistorted near the origin RA/Dec, though the target may appear
#   distorted if it is a very large angular distance from the origin.
ax = body.plot_wireframe_angular(origin_ra=ra0, origin_dec=dec0)
x, y = zip(
    *[
        body.radec2angular(ra, dec, origin_ra=ra0, origin_dec=dec0)
        for ra, dec in zip(ras, decs)
    ]
)
ax.plot(x, y, color='r')
plt.show()

image image image

If this still doesn't exactly do what you want, you may be able to use the values generated by PlanetMapper to save you some effort rather than completely developing everything from scratch. For example, you can use functions like Body.limb_radec to get the coordinates that you can then use as needed. The source for the wireframe plots may also give you some hints for where the different plot elements come from https://github.com/ortk95/planetmapper/blob/main/planetmapper/body.py#L2293-L2420.

I hope this helpful, and do let me know if any of this doesn't make sense, or if you have any other questions!

goatsneez commented 7 months ago

Thanks, this is very useful !

ortk95 commented 7 months ago

Great! I'll close this issue now, but feel free to reopen it/open new issues if you have any more questions. And thanks again for bringing this to my attention!