Create geo-converter.py

geo-converter.py

@@ -0,0 +1,93 @@
+import math
+
+def gcj02_to_wgs84(gcj_lon, gcj_lat):
+    """
+    GCJ-02坐标转换为WGS-84坐标（近似反向转换）
+    """
+    # 迭代逼近法
+    wgs_lon, wgs_lat = gcj_lon, gcj_lat
+    for i in range(30):
+        tmp_lon, tmp_lat = wgs84_to_gcj02(wgs_lon, wgs_lat)
+        wgs_lon = wgs_lon - (tmp_lon - gcj_lon)
+        wgs_lat = wgs_lat - (tmp_lat - gcj_lat)
+    return wgs_lon, wgs_lat
+
+def wgs84_to_gcj02(wgs_lon, wgs_lat):
+    """
+    WGS-84坐标转换为GCJ-02坐标（火星坐标系）
+    
+    Args:
+        wgs_lon (float): WGS-84经度
+        wgs_lat (float): WGS-84纬度
+    
+    Returns:
+        tuple: (gcj_lon, gcj_lat) GCJ-02坐标
+    """
+    
+    # 判断是否在中国范围内
+    def out_of_china(lon, lat):
+        if lon < 72.004 or lon > 137.8347:
+            return True
+        if lat < 0.8293 or lat > 55.8271:
+            return True
+        return False
+    
+    # 如果不在中国范围内，直接返回原坐标
+    if out_of_china(wgs_lon, wgs_lat):
+        return wgs_lon, wgs_lat
+    
+    # 转换参数
+    a = 6378245.0  # 长半轴
+    ee = 0.00669342162296594323  # 扁率
+    
+    def transform_lat(x, y):
+        """纬度转换"""
+        ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(abs(x))
+        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
+        ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0
+        ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0
+        return ret
+    
+    def transform_lon(x, y):
+        """经度转换"""
+        ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(abs(x))
+        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
+        ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0
+        ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0
+        return ret
+    
+    # 计算偏移量
+    d_lat = transform_lat(wgs_lon - 105.0, wgs_lat - 35.0)
+    d_lon = transform_lon(wgs_lon - 105.0, wgs_lat - 35.0)
+    
+    # 纬度弧度
+    rad_lat = wgs_lat / 180.0 * math.pi
+    magic = math.sin(rad_lat)
+    magic = 1 - ee * magic * magic
+    sqrt_magic = math.sqrt(magic)
+    
+    # 调整偏移量
+    d_lat = (d_lat * 180.0) / ((a * (1 - ee)) / (magic * sqrt_magic) * math.pi)
+    d_lon = (d_lon * 180.0) / (a / sqrt_magic * math.cos(rad_lat) * math.pi)
+    
+    # 计算最终坐标
+    gcj_lat = wgs_lat + d_lat
+    gcj_lon = wgs_lon + d_lon
+    
+    return gcj_lon, gcj_lat
+
+
+if __name__ == "__main__":
+    wgs_lon = 116.397428  # 北京天安门经度
+    wgs_lat = 39.90923   # 北京天安门纬度
+
+    print(f"WGS-84坐标: ({wgs_lon}, {wgs_lat})")
+    
+    gcj_lon, gcj_lat = wgs84_to_gcj02(wgs_lon, wgs_lat)
+    
+    print(f"GCJ-02坐标: ({gcj_lon}, {gcj_lat})")
+    print(f"偏移量: ({gcj_lon - wgs_lon}, {gcj_lat - wgs_lat})")
+
+    wgs_lon, wgs_lat = gcj02_to_wgs84(gcj_lon, gcj_lat)
+
+    print(f"WGS-84坐标(转换): ({wgs_lon}, {wgs_lat})")