嘿，我来帮你搞定这两个地球空间计算的问题——先从地球曲率计算说起，再讲3D球体上的区域框怎么算，全程不碰二维投影，纯球面逻辑～

地球本质是个近似椭球，但大多数场景下我们可以简化为球体（用平均半径6371km）。下面给你两个最常用的曲率计算函数，涵盖理论和实际应用场景：

1.1 球面/椭球面曲率半径计算 #

如果需要精确计算地球表面某点的曲率半径（比如地理测绘、高精度导航场景），可以用这个函数——支持球体和椭球两种模型：

import math

def earth_curvature_radius(is_ellipsoid=False, latitude=0.0):
    """
    计算地球表面某点的曲率半径
    :param is_ellipsoid: 是否采用WGS84椭球模型，默认False（球体模型）
    :param latitude: 纬度（弧度），仅椭球模型需要传入
    :return: 曲率半径（单位：米）
    """
    # 地球平均半径（球体模型，所有点曲率半径一致）
    R_avg = 6371000.0
    if not is_ellipsoid:
        return R_avg
    
    # WGS84标准椭球参数
    a = 6378137.0  # 赤道半径（长半轴）
    b = 6356752.3142  # 极半径（短半轴）
    e_sq = (a**2 - b**2) / a**2  # 第一偏心率平方
    
    # 子午圈（南北方向）曲率半径
    R_m = (a * (1 - e_sq)) / math.pow(1 - e_sq * math.sin(latitude)**2, 1.5)
    # 卯酉圈（东西方向）曲率半径
    R_n = a / math.sqrt(1 - e_sq * math.sin(latitude)**2)
    
    # 返回平均曲率半径（取子午圈和卯酉圈的加权平均）
    return (2 * R_n + R_m) / 3

说明 #

• 球体模型下，所有点的曲率半径都是平均半径，适合大多数非高精度场景；
• 椭球模型下，纬度越高，子午圈曲率半径越小，卯酉圈曲率半径越大，这个函数返回的是两者的加权平均，更贴合实际地表的曲率。

1.2 视线遮挡的曲率计算（实用场景） #

这是日常应用中最常用的曲率计算——比如计算距离你d米远的物体，被地球曲率挡住的高度：

def curvature_obstruction_height(distance, R=6371000.0):
    """
    计算距离d处，地球曲率造成的遮挡高度
    :param distance: 水平距离（单位：米）
    :param R: 地球半径（单位：米），默认用平均半径
    :return: 遮挡高度（单位：米）
    """
    # 公式推导：h = R*(sqrt(1 + (d/R)^2) - 1)
    return R * (math.sqrt(1 + (distance / R)**2) - 1)

举个例子 #

如果站在海平面上看10公里外的船，遮挡高度大概是curvature_obstruction_height(10000)≈7.8米——也就是说，船的底部7.8米会被地球曲率挡住，你只能看到船的上部。

你要求不使用二维投影，直接基于地球的3D球体形态计算区域框——这里的区域框其实是球面上由两条纬线（上下边界）和两条经线（左右边界）围成的球面四边形，完全贴合球体表面，没有投影形变。

2.1 核心逻辑 #

1. 把中心点的经纬度、增量值从角度转成弧度（三角函数计算要求）；
2. 计算区域的纬度范围：中心点纬度±Δ纬度/2；
3. 计算区域的经度范围：中心点经度±Δ经度/2，然后把经度归一到[-180°, 180°]（或[-π, π]弧度）区间，避免出现超出范围的数值；
4. 生成四个角点（东北、西北、东南、西南）的经纬度，以及对应的3D笛卡尔坐标（方便后续3D空间计算）。

2.2 实现代码（Python） #

import math

def calculate_spherical_bounding_box(center_lat_deg, center_lon_deg, delta_lat_deg, delta_lon_deg, R=6371000.0):
    """
    计算地球3D球体上的区域框，返回四个角点的经纬度（角度）和笛卡尔坐标
    :param center_lat_deg: 中心点纬度（角度，北纬正、南纬负）
    :param center_lon_deg: 中心点经度（角度，东经正、西经负）
    :param delta_lat_deg: 纬度总增量（角度，比如Δy=10表示上下各扩展5°）
    :param delta_lon_deg: 经度总增量（角度，比如Δx=10表示左右各扩展5°）
    :param R: 地球半径（单位：米），默认用平均半径
    :return: 字典，包含四个角点的详细信息
    """
    # 角度转弧度
    center_lat = math.radians(center_lat_deg)
    center_lon = math.radians(center_lon_deg)
    delta_lat = math.radians(delta_lat_deg)
    delta_lon = math.radians(delta_lon_deg)
    
    # 计算四个角点的弧度坐标
    ne_lat = center_lat + delta_lat / 2
    ne_lon = center_lon + delta_lon / 2
    nw_lat = center_lat + delta_lat / 2
    nw_lon = center_lon - delta_lon / 2
    se_lat = center_lat - delta_lat / 2
    se_lon = center_lon + delta_lon / 2
    sw_lat = center_lat - delta_lat / 2
    sw_lon = center_lon - delta_lon / 2
    
    # 归一化经度到[-π, π]区间（对应角度[-180°, 180°]）
    def normalize_lon(lon_rad):
        while lon_rad > math.pi:
            lon_rad -= 2 * math.pi
        while lon_rad < -math.pi:
            lon_rad += 2 * math.pi
        return lon_rad
    
    ne_lon = normalize_lon(ne_lon)
    nw_lon = normalize_lon(nw_lon)
    se_lon = normalize_lon(se_lon)
    sw_lon = normalize_lon(sw_lon)
    
    # 辅助函数：弧度转角度
    def rad_to_deg(rad):
        return math.degrees(rad)
    
    # 辅助函数：球面坐标转3D笛卡尔坐标
    def spherical_to_cartesian(lat_rad, lon_rad):
        x = R * math.cos(lat_rad) * math.cos(lon_rad)
        y = R * math.cos(lat_rad) * math.sin(lon_rad)
        z = R * math.sin(lat_rad)
        return (round(x, 2), round(y, 2), round(z, 2))
    
    # 组装结果
    bounding_box = {
        "northeast": {
            "lat_deg": round(rad_to_deg(ne_lat), 4),
            "lon_deg": round(rad_to_deg(ne_lon), 4),
            "cartesian": spherical_to_cartesian(ne_lat, ne_lon)
        },
        "northwest": {
            "lat_deg": round(rad_to_deg(nw_lat), 4),
            "lon_deg": round(rad_to_deg(nw_lon), 4),
            "cartesian": spherical_to_cartesian(nw_lat, nw_lon)
        },
        "southeast": {
            "lat_deg": round(rad_to_deg(se_lat), 4),
            "lon_deg": round(rad_to_deg(se_lon), 4),
            "cartesian": spherical_to_cartesian(se_lat, se_lon)
        },
        "southwest": {
            "lat_deg": round(rad_to_deg(sw_lat), 4),
            "lon_deg": round(rad_to_deg(sw_lon), 4),
            "cartesian": spherical_to_cartesian(sw_lat, sw_lon)
        }
    }
    
    return bounding_box

2.3 调用示例 #

比如中心点是北纬30°、东经120°，纬度增量10°，经度增量10°：

bbox = calculate_spherical_bounding_box(30, 120, 10, 10)
print(bbox["northeast"])
# 输出：{'lat_deg': 35.0, 'lon_deg': 125.0, 'cartesian': (3185546.02, 5518143.01, 3337890.04)}

2.4 关键提醒 #

• 为什么不用二维投影？因为任何二维投影都会带来形变（比如墨卡托投影在高纬度会拉伸），而这个方法直接计算球面上的点，得到的区域是真实的球体表面区域，边界是经线和纬线的弧段；
• 经度归一化很重要：经度是循环的（比如190°西经=170°东经），必须把计算后的经度归一到标准区间，避免出现异常数值；
• 笛卡尔坐标的作用：如果需要判断某个点是否在这个区域内、计算区域的球面面积等，用笛卡尔坐标做向量运算会更方便。

内容的提问来源于stack exchange，提问作者four-eyes