Instantly share code, notes, and snippets.

Embed
What would you like to do?
Mercator to BD09, BD09 to GCJ03
#!/usr/bin/python
#-*- coding=utf-8 -*-
"""用于和地理信息有关的数据处理
"""
__author__ = ['"wuyadong" <wuyadong@tigerknows.com>']
import math
MCBAND = (12890594.86, 8362377.87, 5591021, 3481989.83, 1678043.12, 0)
MC2LL = ([1.410526172116255e-8, 0.00000898305509648872, -1.9939833816331,
200.9824383106796, -187.2403703815547, 91.6087516669843, - 23.38765649603339,
2.57121317296198, -0.03801003308653, 17337981.2],
[-7.435856389565537e-9, 0.000008983055097726239, -0.78625201886289,
96.32687599759846, -1.85204757529826, -59.36935905485877, 47.40033549296737,
-16.50741931063887, 2.28786674699375, 10260144.86],
[-3.030883460898826e-8, 0.00000898305509983578, 0.30071316287616,
59.74293618442277, 7.357984074871, -25.38371002664745, 13.45380521110908,
-3.29883767235584, 0.32710905363475, 6856817.37],
[-1.981981304930552e-8, 0.000008983055099779535, 0.03278182852591, 40.31678527705744,
0.65659298677277, -4.44255534477492, 0.85341911805263, 0.12923347998204,
-0.04625736007561, 4482777.06],
[3.09191371068437e-9, 0.000008983055096812155, 0.00006995724062, 23.10934304144901,
-0.00023663490511, -0.6321817810242, -0.00663494467273, 0.03430082397953,
-0.00466043876332, 2555164.4],
[2.890871144776878e-9, 0.000008983055095805407, -3.068298e-8, 7.47137025468032,
-0.00000353937994, -0.02145144861037, -0.00001234426596, 0.00010322952773,
-0.00000323890364, 826088.5])
X_PI = 3.14159265358979324 * 3000.0 / 180.0
class GISError(Exception):
"""GIS Exception
"""
def convert_MCT_2_BD09(lon, lat):
"""将墨卡托坐标转换成BD09
Args:
lon: float, 经度
lat: float, 维度
Returns:
(x, y): tuple, 经过转换的x, y
"""
ax = None
# 获取常量ax
for j in xrange(len(MCBAND)):
if lat >= MCBAND[j]:
ax = MC2LL[j]
break
if ax is None:
raise GISError("error lat:%s" % lat)
e = ax[0] + ax[1] * abs(lon)
i = abs(lat) / ax[9]
aw = ax[2] + ax[3] * i + ax[4] * i * i + ax[5] * i * i * i +\
ax[6] * i * i * i * i + ax[7] * i * i * i * i * i + ax[8] * i * i * i * i * i * i
if lon < 0:
e *= -1
if lat < 0:
aw *= -1
return e, aw
def convert_BD09_2_GCJ03(lon, lat):
"""坐标转换,将BD09转成GCJ03
Args:
lon: float, 经度
lat: float, 维度
Returns:
(x, y): tuple, 转换后的结果
"""
x = lon - 0.0065
y = lat - 0.006
z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * X_PI)
theta = math.atan2(y, x) - 0.000003 * math.cos(x * X_PI)
gg_lon = z * math.cos(theta)
gg_lat = z * math.sin(theta)
return gg_lon, gg_lat
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment