计算地球上空间两点距离

计算地球上空间两点距离

geopy.distance.great_circle和scipy.spatial.cKDTree两种方法：

• geopy.distance.great_circle:

import numpy as np
import pandas as pd
from geopy.distance import great_circle

# Load the data
ams_data = pd.read_csv('cutoff_energy1days_calc.csv')
nm0_data = pd.read_csv('data_v2_check_calc.csv', index_col=0)
nm_data = nm0_data.T.reset_index()

# Function to find the closest point
def find_closest_point(nm_point, ams_data):
    min_distance = float('inf')  # Initialize min_distance with infinity
    closest_index = None  # Initialize closest_index with None
    # Loop through all points in ams_data
    for j in range(len(ams_data)):
        point_a = (nm_point['latitude_degree'], nm_point['longitude_degree'])
        point_b = (ams_data.iloc[j]['latitude_degree'], ams_data.iloc[j]['longitude_degree'])
        distance = great_circle(point_a, point_b).kilometers
        if distance < min_distance:
            min_distance = distance
            closest_index = j
    return closest_index, min_distance

# Iterate over nm_data and find the closest point from ams_data
for i, nm_point in nm_data.iterrows():
    closest_index, min_distance = find_closest_point(nm_point, ams_data)
    print(nm_data['index'][i], ams_data.iloc[closest_index]['utime'], min_distance, ams_data.iloc[closest_index]['CutoffRigidity_GV'])

# Note: You may need to adjust the column names based on the actual column names in your CSV files

• scipy.spatial.cKDTree:

import numpy as np
import pandas as pd
from scipy.spatial import cKDTree

# Load the data
ams_data = pd.read_csv('cutoff_energy20days_calc.csv')
nm0_data = pd.read_csv('data_v2_check_calc.csv', index_col=0)
nm_data = nm0_data.T.reset_index()

# Convert degrees to radians for the spatial tree
def degrees_to_radians(dataframe, lat_col, lon_col):
    dataframe[lat_col] = np.radians(dataframe[lat_col])
    dataframe[lon_col] = np.radians(dataframe[lon_col])
    
degrees_to_radians(ams_data, 'latitude_degree', 'longitude_degree')
degrees_to_radians(nm_data, 'latitude_degree', 'longitude_degree')

# Build the spatial tree
ams_coords = ams_data[['latitude_degree', 'longitude_degree']].values
tree = cKDTree(ams_coords)

# Function to find the closest point using cKDTree
def find_closest_point(tree, nm_point):
    distance, index = tree.query([nm_point['latitude_degree'], nm_point['longitude_degree']])
    return index, distance

# Iterate over nm_data and find the closest point from ams_data using cKDTree
for i, nm_point in nm_data.iterrows():
    closest_index, min_distance = find_closest_point(tree, nm_point)
    # Convert the distance from radians to kilometers assuming Earth's radius of 6371 km
    min_distance_km = min_distance * 6371
    print(nm_data['index'][i], ams_data.iloc[closest_index]['utime'], min_distance_km,ams_data.iloc[closest_index]['CutoffRigidity_GV'])

# Note: Make sure to use radians for latitude and longitude when working with cKDTree