Question about Calculating Tick Coordinates in Astropy

Astropy
, ,
1

Hi all,
I am creating an application in which I need to display HMI magnetograms in CEA projections. Because Matplotlib is not an option for this project, I need to plot everything manually with a different package.

To be able to graph these images, I need to label the ticks on the axes, but I’m having trouble finding the correct tick locations. I need pixel coordinates to graph ticks, but I want ticks separated at intervals based on the world coordinates. For example, if I needed a tick at a latitude of 30 degrees, I can’t just convert any arbitrary coordinate along the 30 degree latitude line to pixel coordinates, because the line is curved in the image. So I need to be able to find the specific pixel coordinates of a point whose latitude is at 30 degrees but also at the 0 pixel x coordinate. Basically, I’d like to be able to find the coordinates of a point using two different coordinate systems for x and y.

I know it is possible to calculate these manually, but as I am using Astropy wcs for everything else coordinate-related, it would be nice to keep everything consistent. Is there any way to do this?

Thanks in advance

2

Hey @a19231 :wave:

I am not sure I fully grasp your exact question, would I be correct in saying that you want to take a line of constant latitude/longitude (in HGC?) and find the pixel coordinates for that line?

The way that WCSAxes (the wcs aware mpl plugin Astropy provides) works is that it locates the ticks based on world coordinates not pixel coordinates, so it finds the range of world coordinates that intersect with the x spine and spaces ticks out along that axis.

Here is an example of how I would manually draw a grid on a Carrington map:

import numpy as np
import matplotlib.pyplot as plt

import astropy.units as u
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
from astropy.utils.data import download_file
from astropy.wcs.utils import wcs_to_celestial_frame

import sunpy.coordinates


filename = download_file(
    'http://jsoc.stanford.edu/data/hmi/synoptic/hmi.Synoptic_Mr.2191.fits', cache=True)

hdu = fits.open(filename)[0]
header, data = hdu.header, hdu.data
header["CUNIT1"] = "deg"
header["CUNIT2"] = "deg"
header["CDELT2"] = 180/np.pi * header["CDELT2"]

wcs = WCS(header)
frame = wcs_to_celestial_frame(wcs)

# Note, for a more non-linear coordinate space you would want to calculate lat and lon separately.
# I am just being lazy and only using two corners rather than all 4
start_hgc = wcs.pixel_to_world(0, 0)
end_hgc = wcs.pixel_to_world(*np.array(data.shape[::-1])-1)

lon_lines = SkyCoord(np.linspace(start_hgc.lon, end_hgc.lon, 10, endpoint=True)[:, None],
                     np.linspace(start_hgc.lat, end_hgc.lat, 10, endpoint=True),
                     unit=u.deg,
                     frame=frame)
lon_pixels = wcs.world_to_pixel(lon_lines)

lat_lines = SkyCoord(np.linspace(start_hgc.lon, end_hgc.lon, 10, endpoint=True),
                     np.linspace(start_hgc.lat, end_hgc.lat, 10, endpoint=True)[:, None],
                     unit=u.deg,
                     frame=frame)
lat_pixels = wcs.world_to_pixel(lat_lines)


fig, ax = plt.subplots()

ax.imshow(data, origin="lower", interpolation="none")
ax.plot(*lon_pixels, "--", color="orange")
ax.plot(*lat_pixels, "--", color="orange")
plt.show()

If you have other questions please let us know :slight_smile:

1 Like
3

Which gives you this:

Figure_1