I’m new to Astropy and would appreciate some help. I’m trying to mask out pixels in images from a space telescope whose lines of sight pass through e.g. the Lunar or Earth discs.
I’ve found this sunpy example which is quite close to what I’m looking for, but needs to work for the Earth and Moon. Any guidance would be appreciated.
What I have:
image data in ndarray
instrument orientation unit vectors (3x3 matrix in GCRS)
instrument position (3x1 matrix in GCRS) and date
instrument FOV
What I want:
boolean array mask whether each pixel’s LOS intersects a sphere of size r centered around a SkyCoord returned by get_body()
Here are two more pieces of information you’ll need:
How the pixels map to spherical coordinates (in the parlance of WCS, the appropriate projection)
Whether the instrument’s orientation vectors in GCRS are already corrected for stellar aberration (i.e., the shift in the apparent direction of an observed object as a result of observer motion).
At high level, here’s how to accomplish what you want for the Earth disc (and the Moon is just the same*).
Use get_body() to get the location of the Earth in GCRS. Ideally you’d specify the accurate observer location and velocity.
Subtract the observer location from the location of the instrument to get the observer-Earth vector.
Take the arcsine of the ratio of the Earth radius and the observer-Earth distance to get the apparent angular size of the Earth.
For every pixel in your data array, calculate the corresponding 2D direction in GCRS for the center of that pixel. This is straightforward if you are working with a constructed WCS. (Of course, adjust your approach if you need something more sophisticated than masking based on pixel centers.)
For every pixel in your data array, calculate the angle between the above 2D direction and the observer-Earth vector (hint: take the dot product)
The mask you want is simply where the above angle array is less than the apparent angular size of the Earth.
* I’ll point out that the position of the Moon can be appreciably inaccurate with Astropy’s built-in ephemeris. I highly recommend you use a JPL ephemeris for Moon calculations. See Solar System Ephemerides — Astropy v5.3.3
Constructing a WCS header means that you can have wcslib compute the pixel->world conversion, and thus avoid re-implementing the transformation code and potentially introducing errors. For the gnomonic projection, a common simplification that people do, sometimes unintentionally, is to use the small-angle approximation, and inaccuracies can creep in.
Besides, having a WCS header means that you can do stuff like use the reproject package to reproject your image into other projections.