step1

当网格号内的link跨越UTM分区的网格时，会遇到如下问题。

问题： 每个UTM分区有独立的中央经线和坐标系，直接转换可能会导致不同区的坐标不一致。简单地使用 utm.from_latlon() 会使得不同UTM分区的坐标不具备直接可比性。

解决方案： 强制转换到一个固定的UTM区。手动选择一个目标UTM区，将所有经纬度坐标都转换为这个区的UTM坐标。

import utm

def from_latlon_fixed(lat, lon, zone_number):
    """
    将经纬度坐标转换到指定的 UTM 区域中。
    lat: 纬度
    lon: 经度
    zone_number: 目标UTM区号
    """
    utm_result = utm.from_latlon(lat, lon, force_zone_number=zone_number)
    return utm_result
    
if __name__ == "__main__":
   start_t = time.time()
   
   fixed_zone_number = 33 
   lati = 41.0 
   longi = -101.0 
   res = from_latlon_fixed(lati, longi, fixed_zone_number) 
   print(f"east = {res[0]}, north = {res[1]}, utm_number = {res[2]}.")      
   
   end_t = time.time()
   print(f"used time: {(end_t - start_t) / 60.0} minutes!")

输出为，

east = -6588362.84346626, north = 18159942.833583914, utm_number = 33.
used time: 0.0 minutes!

step2

计算网格meshid = 605601的中心点经度纬度，

def get_mesh_grid(mesh_id):
    """
    根据标准图幅号算出角点坐标
    return: [min_longi, max_longi, min_lati, max_lati]
    """
    ide = round(float(mesh_id))
    min_x = (int(ide / 100) % 100 * 8 + ide % 10) / 8. + 60
    min_y = (int(ide / 10000) * 8 + int(ide / 10) % 10) / 12.
    max_x = min_x + 1 / 8.
    max_y = round(min_y + 1 / 12., 8)
    min_y = round(min_y, 8)
    return [min_x, min_y, max_x, max_y]

def get_center_of_meshid(meshid):
    """
    根据标准图幅号计算中心点经度纬度
    """
    res = get_mesh_grid(meshid)
    center_longi = (res[0] + res[2]) / 2.0 
    center_lati = (res[1] + res[3]) / 2.0 
    return [center_longi, center_lati]

if __name__ == "__main__":
    meshid = "605601"
    longi, lati = get_center_of_meshid(meshid)
    print(f"longi = {longi}, lati = {lati}.")

输出为，

longi = 116.1875, lati = 40.041666665.
used time: 0.0 minutes!

step3

根据给定的经纬度，计算其对应的UTM分区号是一个相对简单的过程。UTM 系统将地球表面分为60个纵向的分区，每个分区宽度为6°经度。UTM分区号（Zone Number）是基于经度来确定的，公式如下：

UTM Zone Number 的计算公式为：

\text{Zone Number} = \left\lfloor \frac{\text{Longitude} + 180}{6} \right\rfloor + 1

北纬还是南纬会影响 UTM 坐标中的区字母，但不影响区号。

示例 Python 代码

def calculate_utm_zone(longitude, latitude):
    """
    根据经度和纬度计算对应的 UTM 分区号
    """
    # 计算 UTM 分区号
    zone_number = int((longitude + 180) / 6) + 1
    # 北半球和南半球使用不同的区字母
    if latitude >= 0:
        hemisphere = 'Northern Hemisphere'
    else:
        hemisphere = 'Southern Hemisphere'
    return zone_number, hemisphere

if __name__ == "__main__":
    longitude = -74.0060  # 例如：纽约市
    latitude = 40.7128
    zone_number, hemisphere = calculate_utm_zone(longitude, latitude)
    print(f"zone_number = {zone_number}, hemisphere = {hemisphere}.")

运行结果为，

zone_number = 18, hemisphere = Northern Hemisphere.

将一个网格号内的所有link（经度纬度格式）转到同一个平面直角坐标系

step1

step2

step3

示例 Python 代码