高德火星坐标(GCJ-02)转WGS84坐标
1 转换算法
import mathdef gcj02_to_wgs84(lon, lat):"""高德火星坐标(GCJ-02)转WGS84坐标"""a = 6378245.0 ee = 0.00669342162296594323 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.0ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0return retdef 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.0ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0return retdlat = transform_lat(lon - 105.0, lat - 35.0)dlon = transform_lon(lon - 105.0, lat - 35.0)radlat = lat / 180.0 * math.pimagic = math.sin(radlat)magic = 1 - ee * magic * magicsqrtmagic = math.sqrt(magic)dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * math.pi)dlon = (dlon * 180.0) / (a / sqrtmagic * math.cos(radlat) * math.pi)wgs_lat = lat - dlatwgs_lon = lon - dlonreturn wgs_lon, wgs_lat


2 使用Arcpy实现点要素类坐标转换
import arcpy
import mathdef gcj02_to_wgs84(lon, lat):"""高德火星坐标(GCJ-02)转WGS84坐标"""a = 6378245.0 ee = 0.00669342162296594323 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.0ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0return retdef 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.0ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0return retdlat = transform_lat(lon - 105.0, lat - 35.0)dlon = transform_lon(lon - 105.0, lat - 35.0)radlat = lat / 180.0 * math.pimagic = math.sin(radlat)magic = 1 - ee * magic * magicsqrtmagic = math.sqrt(magic)dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * math.pi)dlon = (dlon * 180.0) / (a / sqrtmagic * math.cos(radlat) * math.pi)wgs_lat = lat - dlatwgs_lon = lon - dlonreturn wgs_lon, wgs_lat
input_fc = r"Point"
output_fc = r"Point2"
output_sr = arcpy.SpatialReference(4490) ws = r"D:\admin\Documents\Projects\tj\tj.gdb"arcpy.env.workspace = ws
arcpy.management.CreateFeatureclass(out_path=arcpy.env.workspace,out_name=output_fc,geometry_type="POINT",spatial_reference=output_sr
)
arcpy.management.AddField(output_fc, "NAME", "TEXT")
with arcpy.da.Editor(arcpy.env.workspace) as edit:with arcpy.da.InsertCursor(output_fc, ["SHAPE@", "NAME"]) as insert_cursor:with arcpy.da.SearchCursor(input_fc, ["SHAPE@", "NAME"]) as cursor:for row in cursor:original_point = row[0]x = original_point.firstPoint.Xy = original_point.firstPoint.Yconverted_x, converted_y = gcj02_to_wgs84(x, y)new_point = arcpy.Point(converted_x, converted_y)new_geometry = arcpy.PointGeometry(new_point, output_sr)insert_cursor.insertRow([new_geometry, row[1]])print("坐标转换完成!")
ut_sr)insert_cursor.insertRow([new_geometry, row[1]])print("坐标转换完成!")