优化从 Astropy 的地心真黄道天空坐标中提取纬度/经度

问题描述 投票：0回答：1

我面临着一个挑战，这个挑战应该很简单，但事实证明 Astropy 库相当复杂。具体来说，我需要经常计算已转换为 

SkyCoord
框架的 
GeocentricTrueEcliptic
对象的经度和纬度。这是我的代码的相关部分：

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)
由于不必要的计算与我的用例不符，

访问 

coo.lon
和 
coo.lat
被证明非常耗时。这个问题变得更加严重，因为在我的模块的每个测试阶段都需要重复计算数百万次，从而显着影响性能。

经过调查，我发现

coo.lon
和
coo.lat
可能会很慢，因为
SkyCoord
执行到球坐标的转换。我通过直接访问内部属性发现了一种更快的方法（
coo_icrs._sky_coord_frame._data
）：

%timeit coo_icrs._sky_coord_frame._data._lon
%timeit coo_icrs.ra
# 46.1 ns ± 0.203 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)
# 6.21 µs ± 69.8 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

但是，此快捷方式不适用于 

coo
，因为它的坐标仅以 XYZ 笛卡尔格式提供：

coo._sky_coord_frame._data
# <CartesianRepresentation (x, y, z) [dimensionless]
#    (-0.83586289, 0.53450713, 0.12504145)>

直接提取值（

._x.value
、
._y.value
和
._z.value
）并手动将其转换为经度和纬度仍然比使用
coo.lon
和
coo.lat
更快，但这种手动转换的计算量也很大。

我正在寻找一种方法来更有效地从 

SkyCoord
帧中的 
GeocentricTrueEcliptic
对象检索经度和纬度值，避免额外的计算开销。有没有更快、更直接的方法来获取这些值？

python numpy numba jit astropy
1个回答
2
投票

以下答案是我的博客增强从 Astropy 的地心真黄道框架中提取纬度和经度的性能的简洁摘录。

下面介绍的三个函数使用 

x
对象中的 
y
、
z
和 
astropy.coordinates.sky_coordinate.SkyCoord
计算纬度和经度（以弧度表示）。

有一个 

numpy
选项，以及两个 
numba
 增强版本 - 一个使用 JIT，另一个添加并行处理。

结果总结

此分析表明，使用 Numba 增强函数计算球面坐标时性能有所提高。

  • 详细的时间测量：

    • 直接访问 
      .lat
      和 
      .lon
      属性大约需要 
      17.73 µs
      （
      8.89 µs
      为 
      .lat
      ，
      8.84 µs
      为 
      .lon
      ）。
    • 直接检索 
      x
      、
      y
      和 
      z
      值（分别为 
      467 ns
      、
      450 ns
      和 
      458 ns
      ）总和约为 
      1.375 µs
      。
    • cartesian_to_spherical_numba
      计算花费了
      267 ns
      。

  • 将直接值检索和计算的时间结合起来：

    • Numba 增强计算的总时间 (
      cartesian_to_spherical_numba
      )：
      1.375 µs
      （值检索）+ 
      267 ns
      = 
      1.642 µs
      。
    • 这比 
      .lat
      和 
      .lon
      (
      17.73 µs / 1.642 µs ≈ 10.8
      ) 的直接访问方法快约 10.8 倍。

这些结果证明了使用 JIT 编译所取得的成果。如果使用坐标数组，可以通过并行化进一步增强这一点。

相关进口

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：3.12.0
  • Astropy版本：6.0.1
  • Numba 版本：0.59.1
  • Numpy 版本：1.26.4

初始代码

# 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
各种功能

与使用

._sky_coord_frame._data.x.value

 时的
467 ns
 相比，
访问 
6.14 µs
的速度明显更快，仅需要 
.cartesian.x.value
。但是，在属性或方法的名称中使用单个下划线 (
_
) 前缀表示它们旨在“受保护”或供内部使用。请参阅此问题、私有变量和公共与非公共成员了解更多详细信息。

%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)

计算经度和纬度的代码

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)  # Computation of longitude.
    hypot = np.sqrt(x**2 + y**2)  # Computation of the hypotenuse for latitude calculation.
    lat = np.arctan2(z, hypot)  # 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
    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.
    lat = np.empty_like(x)  # Preallocate the array for latitude.
    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

使用方法

使用 
x
、
y
和 
z
作为 
np.float64

进行测试 
# 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)

使用 
x
、
y
和 
z
作为 
np.ndarray

进行测试

如果 

x
、
y
和 
z
值可以配置为三个数组或坐标，请使用。

# 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)

弧度和度、分、秒

.lat
和 
.lon
返回格式为（度、分、秒）的值，而 numpy 和 numba 函数返回弧度值。以下代码提供了额外的格式选项，时序分析中未考虑这些选项。

将弧度转换为度数

@jit(nopython=True, cache=True)
def cartesian_to_spherical_numba_deg(x, y, z):
    lon = np.arctan2(y, x)
    hypot = np.sqrt(x**2 + y**2)
    lat = np.arctan2(z, hypot)
    return np.degrees(lon), np.degrees(lat)  # return results in degrees

格式化为度、分和秒

# 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')))
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.