gcj02 to wgs84

wwyhkq

于 2024-07-19 22:05:29 发布

阅读量219

点赞数 2

文章标签： java 数据库前端

本文链接：https://blog.csdn.net/wwyhkq/article/details/140560574

版权

1、GCJ02 to WGS84 python

import math

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

def gcj02_to_wgs84(lon, lat):
    ee = 0.00669342162296594323
    a = 6378245.0
    dlat = transform_lat(lon - 105.0, lat - 35.0)
    dlon = transform_lon(lon - 105.0, lat - 35.0)
    radlat = lat / 180.0 * math.pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = 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)
    mglat = lat + dlat
    mglon = lon + dlon
    return lon * 2 - mglon, lat * 2 - mglat

# 示例
gcj02_lon, gcj02_lat = 116.41024449916938, 39.91640428150164
wgs84_lon, wgs84_lat = gcj02_to_wgs84(gcj02_lon, gcj02_lat)
print(f"WGS-84坐标: 经度: {wgs84_lon}, 纬度: {wgs84_lat}")

2、WGS84 TO GCJ02

import math

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

def wgs84_to_gcj02(lon, lat):
    if out_of_china(lon, lat):
        return lon, lat
    ee = 0.00669342162296594323
    a = 6378245.0
    dlat = transform_lat(lon - 105.0, lat - 35.0)
    dlon = transform_lon(lon - 105.0, lat - 35.0)
    radlat = lat / 180.0 * math.pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = 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)
    mglat = lat + dlat
    mglon = lon + dlon
    return mglon, mglat

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

# 示例
wgs84_lon, wgs84_lat = 116.40400191537934, 39.91500120257608
gcj02_lon, gcj02_lat = wgs84_to_gcj02(wgs84_lon, wgs84_lat)
print(f"GCJ-02坐标: 经度: {gcj02_lon}, 纬度: {gcj02_lat}")

3、GCJ02 to WGS84 C++

#include <cmath>
#include <iostream>

// 常量定义
const double PI = 3.14159265358979323846;
const double A = 6378245.0; // 长半轴
const double EE = 0.00669342162296594323; // 扁率

// 判断是否在中国境内
bool outOfChina(double lat, double lon) {
    return (lon < 72.004 || lon > 137.8347) || (lat < 0.8293 || lat > 55.8271);
}

// 转换纬度
double transformLat(double x, double y) {
    double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(fabs(x));
    ret += (20.0 * sin(6.0 * x * PI) + 20.0 * sin(2.0 * x * PI)) * 2.0 / 3.0;
    ret += (20.0 * sin(y * PI) + 40.0 * sin(y / 3.0 * PI)) * 2.0 / 3.0;
    ret += (160.0 * sin(y / 12.0 * PI) + 320 * sin(y * PI / 30.0)) * 2.0 / 3.0;
    return ret;
}

// 转换经度
double transformLon(double x, double y) {
    double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(fabs(x));
    ret += (20.0 * sin(6.0 * x * PI) + 20.0 * sin(2.0 * x * PI)) * 2.0 / 3.0;
    ret += (20.0 * sin(x * PI) + 40.0 * sin(x / 3.0 * PI)) * 2.0 / 3.0;
    ret += (150.0 * sin(x / 12.0 * PI) + 300.0 * sin(x / 30.0 * PI)) * 2.0 / 3.0;
    return ret;
}

// GCJ-02 转 WGS-84
void gcj02ToWgs84(double gcjLat, double gcjLon, double &wgsLat, double &wgsLon) {
    if (outOfChina(gcjLat, gcjLon)) {
        wgsLat = gcjLat;
        wgsLon = gcjLon;
        return;
    }
    double dLat = transformLat(gcjLon - 105.0, gcjLat - 35.0);
    double dLon = transformLon(gcjLon - 105.0, gcjLat - 35.0);
    double radLat = gcjLat / 180.0 * PI;
    double magic = sin(radLat);
    magic = 1 - EE * magic * magic;
    double sqrtMagic = sqrt(magic);
    dLat = (dLat * 180.0) / ((A * (1 - EE)) / (magic * sqrtMagic) * PI);
    dLon = (dLon * 180.0) / (A / sqrtMagic * cos(radLat) * PI);
    wgsLat = gcjLat - dLat;
    wgsLon = gcjLon - dLon;
}

// 示例使用
int main() {
    double gcjLat = 31.2304;
    double gcjLon = 121.4737;
    double wgsLat, wgsLon;
    gcj02ToWgs84(gcjLat, gcjLon, wgsLat, wgsLon);
    std::cout << "GCJ-02: " << std::to_string(gcjLat) << ", " << std::to_string(gcjLon) << std::endl;
    std::cout << "WGS-84: " << std::to_string(wgsLat) << ", " << std::to_string(wgsLon) << std::endl;
    return 0;
}

4、WGS84 TO GCJ02 C++

#include <cmath>

const double PI = 3.14159265358979323846;
const double A = 6378245.0;
const double EE = 0.00669342162296594323;

bool outOfChina(double lat, double lon) {
    // 纬度和经度都在中国边界外
    return (lon < 72.004 || lon > 137.8347 || lat < 0.8293 || lat > 55.8271);
}

double transformLat(double x, double y) {
    double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(fabs(x));
    ret += (20.0 * sin(6.0 * x * PI) + 20.0 * sin(2.0 * x * PI)) * 2.0 / 3.0;
    ret += (20.0 * sin(y * PI) + 40.0 * sin(y / 3.0 * PI)) * 2.0 / 3.0;
    ret += (160.0 * sin(y / 12.0 * PI) + 320 * sin(y * PI / 30.0)) * 2.0 / 3.0;
    return ret;
}

double transformLon(double x, double y) {
    double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(fabs(x));
    ret += (20.0 * sin(6.0 * x * PI) + 20.0 * sin(2.0 * x * PI)) * 2.0 / 3.0;
    ret += (20.0 * sin(x * PI) + 40.0 * sin(x / 3.0 * PI)) * 2.0 / 3.0;
    ret += (150.0 * sin(x / 12.0 * PI) + 300.0 * sin(x / 30.0 * PI)) * 2.0 / 3.0;
    return ret;
}

void wgs84ToGcj02(double wgLat, double wgLon, double &mgLat, double &mgLon) {
    if (outOfChina(wgLat, wgLon)) {
        mgLat = wgLat;
        mgLon = wgLon;
        return;
    }
    double dLat = transformLat(wgLon - 105.0, wgLat - 35.0);
    double dLon = transformLon(wgLon - 105.0, wgLat - 35.0);
    double radLat = wgLat / 180.0 * PI;
    double magic = sin(radLat);
    magic = 1 - EE * magic * magic;
    double sqrtMagic = sqrt(magic);
    dLat = (dLat * 180.0) / ((A * (1 - EE)) / (magic * sqrtMagic) * PI);
    dLon = (dLon * 180.0) / (A / sqrtMagic * cos(radLat) * PI);
    mgLat = wgLat + dLat;
    mgLon = wgLon + dLon;
}

// 示例
#include <iostream>
int main() {
    double wgLat = 39.915, wgLon = 116.404;
    double mgLat, mgLon;
    wgs84ToGcj02(wgLat, wgLon, mgLat, mgLon);
    std::cout << "GCJ-02坐标: 经度: " << mgLon << ", 纬度: " << mgLat << std::endl;
    return 0;
}