Post

Enhancing Performance of Latitude and Longitude Extraction from Astropy's GeocentricTrueEcliptic Frame

Posted Updated
By Trenton McKinney 9 min read

Introduction

The post explores an efficient method to extract longitude and latitude values from an Astropy SkyCoord object, specifically when transformed to the GeocentricTrueEcliptic frame. The standard access method, .lon and .lat, incurs unnecessary computations, slowing down the process. Direct access to internal attributes offers a faster alternative, but it doesn’t apply to transformed coordinates in the GeocentricTrueEcliptic frame. The objective is to develop a method that calculates longitude and latitude values from this frame without the computational overhead, as this operation is performed millions of times in each test run.

Associated question on Stack Overflow: Optimizing Lat/Lon Extraction from Astropy’s GeocentricTrueEcliptic SkyCoord

GitHub Repository

1
2
3
4
5
6
7

from astropy.coordinates import SkyCoord, GeocentricTrueEcliptic, ICRS
import astropy.units as u
from astropy.time import Time

coo_icrs = SkyCoord(150*u.deg, 19*u.deg, frame=ICRS(), obstime=Time(2460791., format="jd"))
coo = coo_icrs.geocentrictrueecliptic
some_function(coo.lon, coo.lat)

Solution

The following code demonstrates methods for converting Cartesian coordinates to spherical coordinates, specifically latitude and longitude, from an Astropy SkyCoord object transformed to the GeocentricTrueEcliptic frame. Three functions are presented: a basic NumPy implementation, and two Numba-enhanced versions—one using the JIT compiler and another adding parallel processing capabilities.

  • @jit: The Numba JIT decorator optimizes Python functions by compiling them into machine code at runtime, significantly improving performance for numerical computations.
  • cache: This attribute caches the compiled function on disk, eliminating the need for recompilation in subsequent executions and reducing startup time.
  • parallel: By enabling this attribute, the function supports automatic parallelization, effectively utilizing multiple CPU cores to enhance execution speed.

Results Summary

This analysis demonstrates enhanced performance when using Numba-enhanced functions for calculating spherical coordinates from an Astropy SkyCoord object in the GeocentricTrueEcliptic frame compared to directly accessing .lat and .lon properties.

  • Detailed Time Measurements:
    • Direct access to .lat and .lon properties took about 17.73 µs (8.89 µs for .lat and 8.84 µs for .lon).
    • Retrieving x, y, and z values directly (467 ns, 450 ns, and 458 ns respectively) summed to approximately 1.375 µs.
    • The cartesian_to_spherical_numba computation took 267 ns.
  • Combining the time for direct value retrieval and computation:
    • Total time for Numba-enhanced computation (cartesian_to_spherical_numba): 1.375 µs (value retrieval) + 267 ns = 1.642 µs.
    • This is about 10.8 times faster than the direct access method for .lat and .lon (17.73 µs / 1.642 µs ≈ 10.8).

These results highlight the significant efficiency gains achieved by using JIT compilation through Numba, especially when combined with direct attribute access to perform necessary computations. The parallelization feature in Numba further enhances this efficiency, particularly beneficial for large datasets or operations requiring high computational resources.

Relevant Imports

1
2
3
4
5
6

from astropy.coordinates import SkyCoord, GeocentricTrueEcliptic, ICRS, Angle
import astropy.units as u
from astropy.time import Time
import numpy as np
from numba import jit, prange
from typing import Tuple, Union
  • Python Version: 3.12.3
  • Astropy Version: 6.0.1
  • Numba Version: 0.59.1
  • Numpy Version: 1.26.4

Initial Code

1
2
3
4
5
6

# International Celestial Reference System coordinates
coo_icrs = SkyCoord(150*u.deg, 19*u.deg, frame=ICRS(), obstime=Time(2460791., format="jd"))
# Ecliptic coordinates
ecliptic_coords = coo_icrs.transform_to(GeocentricTrueEcliptic)

print(ecliptic_coords.lat, ecliptic_coords.lon)

[6d21m12.19206726s] [145d28m31.99626325s]

timeit for various functions

Accessing ._sky_coord_frame._data.x.value is significantly faster, taking only 467 ns, compared to 6.14 µs when using .cartesian.x.value. However, in Python, the use of a single underscore (_) prefix in the names of attributes or methods typically denotes that they are intended to be “protected” or for internal use. This naming convention is a common practice among Python developers to indicate that these attributes and methods are not part of the public API of a class or module and should be used with caution.

  • Stability and Compatibility: Using protected members can reduce software stability and maintainability. If the Astropy library refactors or removes these attributes, dependent external code may break.
  • Documentation and Support: Protected members are not usually part of the public API and are often undocumented. Understanding their use and implications typically requires examining the source code or consulting the community.
  • Best Practices: It’s advisable to use public methods and properties provided by the library, which are supported and expected to remain stable. If essential functionality is only available via a protected member, consider contacting the maintainers for a public solution or alternative recommendations.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

%timeit ecliptic_coords.lat
8.89 µs ± 132 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
%timeit ecliptic_coords.lon
8.84 µs ± 413 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
%timeit ecliptic_coords._sky_coord_frame._data
76.5 ns ± 0.507 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
%timeit ecliptic_coords._sky_coord_frame._data.xyz.value
15.8 µs ± 351 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
%timeit ecliptic_coords._sky_coord_frame._data.x.value
467 ns ± 10.4 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
%timeit ecliptic_coords._sky_coord_frame._data.y.value
450 ns ± 3.89 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
%timeit ecliptic_coords._sky_coord_frame._data.z.value
458 ns ± 6.94 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
%timeit ecliptic_coords.cartesian.xyz.value
25.1 µs ± 409 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
%timeit ecliptic_coords.cartesian.x.value
6.14 µs ± 215 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Code to Calculate Longitude and Latitude

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

def cartesian_to_spherical_numpy(
    x: Union[np.ndarray, np.float64],
    y: Union[np.ndarray, np.float64],
    z: Union[np.ndarray, np.float64]
    ) -> Tuple[np.ndarray, np.ndarray]:
    """
    Convert Cartesian coordinates to spherical coordinates using NumPy.

    Args:
    x (Union[np.ndarray, np.float64]): X-coordinates in Cartesian.
    y (Union[np.ndarray, np.float64]): Y-coordinates in Cartesian.
    z (Union[np.ndarray, np.float64]): Z-coordinates in Cartesian.

    Returns:
    Tuple[np.ndarray, np.ndarray]: Tuple of longitude and latitude arrays in radians.
    """
    lon = np.arctan2(y, x)  # Compute longitude using arctan2 for proper handling of quadrant.
    hypot = np.sqrt(x**2 + y**2)  # Compute the hypotenuse of x and y to use for latitude calculation.
    lat = np.arctan2(z, hypot)  # Compute latitude as the angle from the positive z-axis.
    return lon, lat


@jit(nopython=True, cache=True)
def cartesian_to_spherical_numba(
    x: Union[np.ndarray, np.float64],
    y: Union[np.ndarray, np.float64],
    z: Union[np.ndarray, np.float64]
    ) -> Tuple[np.ndarray, np.ndarray]:
    """
    Convert Cartesian coordinates to spherical coordinates using Numba for JIT compilation to speed up the computation.

    Args:
    x (Union[np.ndarray, np.float64]): X-coordinates in Cartesian.
    y (Union[np.ndarray, np.float64]): Y-coordinates in Cartesian.
    z (Union[np.ndarray, np.float64]): Z-coordinates in Cartesian.

    Returns:
    Tuple[np.ndarray, np.ndarray]: Tuple of longitude and latitude arrays in radians.
    """
    lon = np.arctan2(y, x)  # Numba-enhanced computation of longitude.
    hypot = np.sqrt(x**2 + y**2)  # Numba-enhanced computation of the hypotenuse for latitude calculation.
    lat = np.arctan2(z, hypot)  # Numba-enhanced computation of latitude.
    return lon, lat


@jit(nopython=True, parallel=True, cache=True)
def cartesian_to_spherical_numba_parallel(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """
    Convert Cartesian coordinates to spherical coordinates with both parallel processing and JIT compilation using Numba.
    This version leverages multiple CPU cores to handle large arrays more efficiently.

    Args:
    x (np.ndarray): X-coordinates in Cartesian.
    y (np.ndarray): Y-coordinates in Cartesian.
    z (np.ndarray): Z-coordinates in Cartesian.

    Returns:
    Tuple[np.ndarray, np.ndarray]: Tuple of longitude and latitude arrays in radians.
    """
    lon = np.empty_like(x)  # Preallocate the array for longitude to enhance performance.
    lat = np.empty_like(x)  # Preallocate the array for latitude to enhance performance.
    for i in prange(x.size):  # Use parallel range (prange) to distribute loop iterations across multiple threads.
        lon[i] = np.arctan2(y[i], x[i])  # Compute longitude for each element in parallel.
        hypot = np.sqrt(x[i]**2 + y[i]**2)  # Compute the hypotenuse for each element in parallel.
        lat[i] = np.arctan2(z[i], hypot)  # Compute latitude for each element in parallel.
    return lon, lat

Usage

Testing with x, y, and z as np.float64

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# Ensure Numba function is compiled before timing - only used for testing - not required when using functions in a loop
_ = cartesian_to_spherical_numba(np.array([1.0]), np.array([1.0]), np.array([1.0]))

# extract x, y, and z values as numpy.float64
x_ecl = ecliptic_coords._sky_coord_frame._data.x.value
y_ecl = ecliptic_coords._sky_coord_frame._data.y.value
z_ecl = ecliptic_coords._sky_coord_frame._data.z.value

# time the non-parallel functions
%timeit cartesian_to_spherical_numpy(x_ecl, y_ecl, z_ecl)
%timeit cartesian_to_spherical_numba(x_ecl, y_ecl, z_ecl)

2.5 µs ± 749 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
267 ns ± 5.27 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

Testing with x, y, and z as np.ndarray

Use if x, y, and z values can be configured into three arrays of length n.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# Generate test data
n = 1000000
x = np.random.uniform(-1, 1, n)
y = np.random.uniform(-1, 1, n)
z = np.random.uniform(-1, 1, n)

# Ensure Numba function is compiled before timing
_ = cartesian_to_spherical_numba(np.array([1.0]), np.array([1.0]), np.array([1.0]))
_ = cartesian_to_spherical_numba_parallel(np.array([1.0]), np.array([1.0]), np.array([1.0]))

# time all the functions
%timeit cartesian_to_spherical_numpy(x, y, z)
%timeit cartesian_to_spherical_numba(x, y, z)
%timeit cartesian_to_spherical_numba_parallel(x, y, z)

43.3 ms ± 562 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
38 ms ± 302 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
4.99 ms ± 194 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Radians and Degree, Minute, Seconds

In astronomical data handling, particularly with Astropy, angles are often displayed in degrees, minutes, and seconds. However, standard mathematical functions, such as those implemented using Numba, return angles in radians. To align these formats and make Numba function results comparable with Astropy’s coordinate system, it’s necessary to convert radians into degrees. Further conversion into a more readable format may also be required. These conversion options are not included in the timing analysis.

Converting Radians to Degrees

1
2
3
4
5
6
7
8
9

@jit(nopython=True, parallel=True)
def cartesian_to_spherical_numba_parallel(x, y, z):
    lon = np.empty_like(x)
    lat = np.empty_like(x)
    for i in prange(x.size):
        lon[i] = np.arctan2(y[i], x[i])
        hypot = np.sqrt(x[i]**2 + y[i]**2)
        lat[i] = np.arctan2(z[i], hypot)
    return np.degrees(lon), np.degrees(lat)

Formatting as Degrees, Minutes, and Seconds

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# Example radians-to-degree conversion
lon_rad = np.array([2.61799388])  # longitude in radians
lat_rad = np.array([0.33161256])  # latitude in radians

# Convert radians to degrees
lon_deg, lat_deg = np.degrees(lon_rad), np.degrees(lat_rad)

# Create Angle objects and format them
lon_angle = Angle(lon_deg, u.degree)
lat_angle = Angle(lat_deg, u.degree)

# Print formatted angles
print(lon_angle.to_string(unit=u.degree, sep=('d', 'm', 's')))
print(lat_angle.to_string(unit=u.degree, sep=('d', 'm', 's')))

References

Astropy, Python, Performance Optimization
python astropy numba latitude longitude performance optimization
This post is licensed under CC BY 4.0 by the author.