在 utils 子包中增加经纬度计算方位角、距离的工具函数，方便直接读取产品后插值到指定经纬度

cinrad/_utils.pyx

@@ -2,12 +2,13 @@
 cimport numpy as np
 import numpy as np
 cimport cython
-from libc.math cimport sin
+from libc.math cimport sin, cos, atan2, sqrt
 
-cdef double deg2rad, vil_const
+cdef double M_PI, deg2rad, vil_const
 cdef int rm = 8500
 
-deg2rad = 3.141592653589793 / 180
+M_PI = 3.141592653589793
+deg2rad =  M_PI / 180
 vil_const = 3.44e-6
 
 cdef np.ndarray height(np.ndarray distance, double elevation, double radarheight):
@@ -113,4 +114,42 @@ def echo_top(double[:, :, ::1] ref, double[:, ::1] distance, double[::1] elev,
                     w1 = (z1 - threshold) / (z1 - z2)
                     w2 = 1 - w1
                     et[i][j] = w1 * h2 + w2 * h1
-    return np.asarray(et)
+    return np.asarray(et)
+
+@cython.boundscheck(False)
+@cython.cdivision(True)
+@cython.wraparound(False)
+def calculate_distance(const double lon1_rad, const double lon2_rad, const double lat1_rad, const double lat2_rad):
+    cdef double R = 6371000.0
+    cdef double dlat = lat2_rad - lat1_rad
+    cdef double dlon = lon2_rad - lon1_rad
+    cdef double a = (sin(dlat / 2) * sin(dlat / 2) +  
+                    cos(lat1_rad) * cos(lat2_rad) *
+                    sin(dlon / 2) * sin(dlon / 2))
+    cdef double c = 2 * atan2(sqrt(a), sqrt(1 - a)) 
+    cdef double distance = R * c
+    return distance
+
+@cython.boundscheck(False)
+@cython.cdivision(True)
+@cython.wraparound(False)
+def calculate_azimuth(const double lon1_rad, const double lon2_rad, const double lat1_rad, const double lat2_rad):
+    cdef double dlon = lon2_rad - lon1_rad
+    cdef double dx = cos(lat2_rad) * sin(dlon)
+    cdef double dy = (cos(lat1_rad) * sin(lat2_rad) -
+                     sin(lat1_rad) * cos(lat2_rad) * cos(dlon))
+    cdef double azimuth = atan2(dx, dy) 
+    azimuth = (azimuth if azimuth >= 0 else azimuth + 2 * M_PI) 
+    return azimuth
+
+@cython.boundscheck(False)
+@cython.cdivision(True)
+@cython.wraparound(False)
+def geo_to_polar(double lon, double lat, double center_lon, double center_lat):
+    cdef double lon_rad = lon*deg2rad
+    cdef double lat_rad = lat*deg2rad
+    cdef double center_lon_rad = center_lon*deg2rad
+    cdef double center_lat_rad = center_lat*deg2rad
+    cdef double azimuth = calculate_azimuth(center_lon_rad, lon_rad, center_lat_rad, lat_rad)
+    cdef double distance = calculate_distance(center_lon_rad, lon_rad, center_lat_rad, lat_rad) / 1e3
+    return azimuth, distance

cinrad/utils.py

@@ -4,6 +4,7 @@
 from typing import Union, Any
 
 import numpy as np
+import math
 
 from cinrad.constants import deg2rad, vil_const
 from cinrad.projection import height
@@ -146,9 +147,90 @@ def echo_top_py(
                     et[i][j] = w1 * h2 + w2 * h1
     return et
 
+def calculate_distance_py(lon1, lat1, lon2, lat2):
+    """
+    Calculate the great-circle distance between two points on the Earth's surface.
+
+    This function uses the Haversine formula which takes into account the spherical shape of the Earth
+    approximated as a sphere with an average radius.
+
+    Args:
+        lon1 (float): Longitude of the first point.
+        lat1 (float): Latitude of the first point.
+        lon2 (float): Longitude of the second point.
+        lat2 (float): Latitude of the second point.
+
+    Returns:
+        float: The distance between the two points in meters.
+    """
+    R = 6371000  # Earth's average radius in meters
+    # Convert latitude and longitude values from degrees to radians to be used in trigonometric calculations
+    lat1_rad = math.radians(lat1)
+    lat2_rad = math.radians(lat2)
+    lon1_rad = math.radians(lon1)
+    lon2_rad = math.radians(lon2)
+    dlat = lat2_rad - lat1_rad
+    dlon = lon2_rad - lon1_rad
+    # Calculate intermediate value 'a' as per the Haversine formula
+    a = (math.sin(dlat / 2) ** 2 +
+         math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(dlon / 2) ** 2)
+    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
+    distance = R * c
+    return distance
+
+def calculate_azimuth_py(lon1, lat1, lon2, lat2):
+    """
+    Calculate the azimuth (bearing) from the first point to the second point on the Earth's surface.
+
+    The azimuth is the angle measured clockwise from the north direction.
+
+    Args:
+        lon1 (float): Longitude of the starting point.
+        lat1 (float): Latitude of the starting point.
+        lon2 (float): Longitude of the destination point.
+        lat2 (float): Latitude of the destination point.
+
+    Returns:
+        float: The azimuth in radians.
+    """
+    lon1_rad = math.radians(lon1)
+    lon2_rad = math.radians(lon2)
+    lat1_rad = math.radians(lat1)
+    lat2_rad = math.radians(lat2)
+    # Calculate the differences in x and y directions on the spherical coordinate system
+    dx = math.cos(lat2_rad) * math.sin(lon2_rad - lon1_rad)
+    dy = (math.cos(lat1_rad) * math.sin(lat2_rad) -
+          math.sin(lat1_rad) * math.cos(lat2_rad) * math.cos(lon2_rad - lon1_rad))
+    azimuth = math.atan2(dx, dy)
+    azimuth = math.degrees(azimuth)
+    # Adjust the azimuth to be in the range of 0 to 360 degrees
+    azimuth %= 360  
+    return math.radians(azimuth)
+
+def geo_to_polar_py(lon, lat, center_lon, center_lat):
+    """
+    Convert geographical coordinates (longitude and latitude) to polar coordinates (azimuth and distance)
+    relative to a given center point.
+
+    Args:
+        lon (float): Longitude of the point to be converted.
+        lat (float): Latitude of the point to be converted.
+        center_lon (float): Longitude of the center reference point.
+        center_lat (float): Latitude of the center reference point.
+
+    Returns:
+        tuple: A tuple containing the azimuth (in radians) and the distance (in kilometers).
+    """
+    azimuth = calculate_azimuth_py(center_lon, center_lat, lon, lat)
+    distance = calculate_distance_py(center_lon, center_lat, lon, lat) / 1e3  # Convert distance to kilometers
+    return azimuth, distance
+
 try:
     from cinrad._utils import *
 except ImportError:
     # When the C-extension doesn't exist, define the functions in Python.
     echo_top = echo_top_py
-    vert_integrated_liquid = vert_integrated_liquid_py
+    vert_integrated_liquid = vert_integrated_liquid_py
+    calculate_distance = calculate_distance_py
+    calculate_azimuth = calculate_azimuth_py
+    geo_to_polar = geo_to_polar_py