Python 实现常见的坐标系之间的转换

常见的一些坐标系之间的转换
话不多说,直接上代码

import mathx_pi = 3.14159265358979324 * 3000.0 / 180.0pi = 3.1415926535897932384626  # πa = 6378245.0  # 长半轴ee = 0.00669342162296594323  # 偏心率平方class CoordinateTransform:    @classmethod    def gcj02_to_bd09(cls, lng, lat): """ 火星坐标系(GCJ-02)转百度坐标系(BD-09) 谷歌、高德——>百度 :param lng:火星坐标经度 :param lat:火星坐标纬度 :return: """ z = math.sqrt(lng * lng + lat * lat) + 0.00002 * math.sin(lat * x_pi) theta = math.atan2(lat, lng) + 0.000003 * math.cos(lng * x_pi) bd_lng = z * math.cos(theta) + 0.0065 bd_lat = z * math.sin(theta) + 0.006 return [bd_lng, bd_lat]    @classmethod    def bd09_to_gcj02(cls, bd_lng, bd_lat): """ 百度坐标系(BD-09)转火星坐标系(GCJ-02) 百度——>谷歌、高德 :param bd_lat:百度坐标纬度 :param bd_lon:百度坐标经度 :return:转换后的坐标列表形式 """ x = bd_lng - 0.0065 y = bd_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) gcj_lng = z * math.cos(theta) gcj_lat = z * math.sin(theta) return [gcj_lng, gcj_lat]    @classmethod    def wgs84_to_gcj02(cls, lng, lat): """ WGS84转GCJ02(火星坐标系) 天地图> 高德谷歌 :param lng:WGS84坐标系的经度 :param lat:WGS84坐标系的纬度 :return: """ if cls.out_of_china(lng, lat):  # 判断是否在国内     return [lng, lat] dlat = cls.transform_lat(lng - 105.0, lat - 35.0) dlng = cls.transform_lng(lng - 105.0, lat - 35.0) radlat = lat / 180.0 * pi magic = math.sin(radlat) magic = 1 - ee * magic * magic sqrtmagic = math.sqrt(magic) dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi) dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi) mglat = lat + dlat mglng = lng + dlng return [mglng, mglat]    @classmethod    def gcj02_to_wgs84(cls, lng, lat): """ GCJ02(火星坐标系)转GPS84 :param lng:火星坐标系的经度 :param lat:火星坐标系纬度 :return: """ if cls.out_of_china(lng, lat):     return [lng, lat] dlat = cls.transform_lat(lng - 105.0, lat - 35.0) dlng = cls.transform_lng(lng - 105.0, lat - 35.0) radlat = lat / 180.0 * pi magic = math.sin(radlat) magic = 1 - ee * magic * magic sqrtmagic = math.sqrt(magic) dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi) dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi) mglat = lat + dlat mglng = lng + dlng return [lng * 2 - mglng, lat * 2 - mglat]    @classmethod    def bd09_to_wgs84(cls, bd_lng, bd_lat): """ 百度09坐标系转 GPS84坐标系 :param bd_lon: :param bd_lat: :return: """ lon, lat = cls.bd09_to_gcj02(bd_lng, bd_lat) return cls.gcj02_to_wgs84(lon, lat)    @classmethod    def wgs84_to_bd09(cls, lng, lat): """ GPS84坐标系转百度09坐标系 :param lon: :param lat: :return: """ lon, lat = cls.wgs84_to_gcj02(lng, lat) return cls.gcj02_to_bd09(lon, lat)    @staticmethod    def transform_lat(lng, lat): ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng)) ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *  math.sin(2.0 * lng * pi)) * 2.0 / 3.0 ret += (20.0 * math.sin(lat * pi) + 40.0 *  math.sin(lat / 3.0 * pi)) * 2.0 / 3.0 ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *  math.sin(lat * pi / 30.0)) * 2.0 / 3.0 return ret    @staticmethod    def transform_lng(lng, lat): ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng)) ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *  math.sin(2.0 * lng * pi)) * 2.0 / 3.0 ret += (20.0 * math.sin(lng * pi) + 40.0 *  math.sin(lng / 3.0 * pi)) * 2.0 / 3.0 ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *  math.sin(lng / 30.0 * pi)) * 2.0 / 3.0 return ret    @classmethod    def out_of_china(cls, lng, lat): """ 判断是否在国内，不在国内不做偏移 :param lng: :param lat: :return: """ return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)    # gcj02经纬度坐标>墨卡托坐标    @classmethod    def lonLat2Mercator(cls, lng: float, lat: float): """ :param longitude: :param latitude: :return: 国测局02坐标(火星坐标) 转 墨卡托坐标x,y """ x = lng * 20037508.34 / 180 y = math.log(math.tan((90 + lat) * math.pi / 360)) / (math.pi / 180) y = y * 20037508.34 / 180 return [x, y]    # 墨卡托坐标>gcj02经纬度坐标    @classmethod    def Mercator2lonLat(cls, mercator_x, mercator_y): x = mercator_x / 20037508.34 * 180 y = mercator_y / 20037508.34 * 180 y = 180 / math.pi * (2 * math.atan(math.exp(y * math.pi / 180)) - math.pi / 2) lng, lat = x, y return [lng, lat]    # 墨卡托坐标>计算瓦片坐标(瓦片编号坐标) y代表地图缩放级别    @classmethod    def wmc2tile(cls, x, y, z=15): """ :param x: mercator_x 坐标 :param y: mercator_y 坐标 :param z: 瓦片坐标缩放层级 :return: """ unit = 40075016.68 / math.pow(2, z) tx = int(((x + 20037508.34) / unit)) ty = int(((20037508.34 - y) / unit)) return [tx, ty, z]if __name__ == '__main__':    cdt = CoordinateTransform()    lng_lat = input(f"请输入经纬度坐标:").strip().strip('\n').split(',')    lng, lat = lng_lat[0], lng_lat[-1]    print(f"原始经纬度坐标:{lng},{lat}")    print(f"转为bd09坐标系:{cdt.gcj02_to_bd09(float(lng), float(lat))}")    print(f"转为wgs84坐标系:{cdt.gcj02_to_wgs84(float(lng), float(lat))}")    mer_coord = cdt.lonLat2Mercator(float(lng), float(lat))    print(f"转为墨卡托坐标:{mer_coord}")    print(f"高德瓦片坐标:{cdt.wmc2tile(mer_coord[0], mer_coord[-1])}")

希望对各位giser有一定帮助