From 246aba1a0def8b1811ef166a63027b8f1a972629 Mon Sep 17 00:00:00 2001 From: Longtsing Date: Thu, 5 Dec 2024 23:50:37 +0800 Subject: [PATCH 1/2] =?UTF-8?q?=E5=9C=A8=20utils=20=E5=AD=90=E5=8C=85?= =?UTF-8?q?=E4=B8=AD=E5=A2=9E=E5=8A=A0=E7=BB=8F=E7=BA=AC=E5=BA=A6=E8=AE=A1?= =?UTF-8?q?=E7=AE=97=E6=96=B9=E4=BD=8D=E8=A7=92=E3=80=81=E8=B7=9D=E7=A6=BB?= =?UTF-8?q?=E7=9A=84=E5=B7=A5=E5=85=B7=E5=87=BD=E6=95=B0=EF=BC=8C=E6=96=B9?= =?UTF-8?q?=E4=BE=BF=E7=9B=B4=E6=8E=A5=E8=AF=BB=E5=8F=96=E4=BA=A7=E5=93=81?= =?UTF-8?q?=E5=90=8E=E6=8F=92=E5=80=BC=E5=88=B0=E6=8C=87=E5=AE=9A=E7=BB=8F?= =?UTF-8?q?=E7=BA=AC=E5=BA=A6?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- cinrad/_utils.pyx | 45 +++++++++++++++++++++++-- cinrad/utils.py | 83 ++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 124 insertions(+), 4 deletions(-) diff --git a/cinrad/_utils.pyx b/cinrad/_utils.pyx index 8103f76..090cf46 100644 --- a/cinrad/_utils.pyx +++ b/cinrad/_utils.pyx @@ -4,10 +4,11 @@ import numpy as np cimport cython from libc.math cimport sin -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) \ No newline at end of file + 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 \ No newline at end of file diff --git a/cinrad/utils.py b/cinrad/utils.py index 6950392..4a06e84 100644 --- a/cinrad/utils.py +++ b/cinrad/utils.py @@ -146,9 +146,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 \ No newline at end of file + vert_integrated_liquid = vert_integrated_liquid_py + calculate_distance = calculate_distance_py + calculate_azimuth = calculate_azimuth_py + geo_to_polar = geo_to_polar_py \ No newline at end of file From 06a6793a9f19a7f6d95814a4b25849189c683d74 Mon Sep 17 00:00:00 2001 From: Longtsing Date: Thu, 5 Dec 2024 23:57:12 +0800 Subject: [PATCH 2/2] =?UTF-8?q?=E5=A2=9E=E5=8A=A0=E4=BA=86=E7=9B=B8?= =?UTF-8?q?=E5=85=B3=E5=BA=93=E7=9A=84=E5=AF=BC=E5=85=A5=E4=BB=A3=E7=A0=81?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- cinrad/_utils.pyx | 2 +- cinrad/utils.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/cinrad/_utils.pyx b/cinrad/_utils.pyx index 090cf46..40e399d 100644 --- a/cinrad/_utils.pyx +++ b/cinrad/_utils.pyx @@ -2,7 +2,7 @@ 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 M_PI, deg2rad, vil_const cdef int rm = 8500 diff --git a/cinrad/utils.py b/cinrad/utils.py index 4a06e84..4aba200 100644 --- a/cinrad/utils.py +++ b/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