地图坐标系基本知识及常用计算

20年笔记迁移

我们常说的坐标系有哪些？

WGS84：为一种大地坐标系，也是目前广泛使用的GPS全球卫星定位系统使用的坐标系。

GCJ02：又称火星坐标系，是由中国国家测绘局制定的地理坐标系统，是由WGS84加密后得到的坐标系。

BD09：为百度坐标系，在GCJ02坐标系基础上再次加密。其中bd09ll表示百度经纬度坐标，bd09mc表示百度墨卡托米制坐标

百度地图使用的是BD09坐标系，若使用非BD09坐标、未经过坐标转换（非BD09转成BD09）直接叠加在地图上，地图展示位置会偏移，因此通过其他坐标（WGS84、GCJ02）调用服务时，需先将其他坐标转换为BD09。

百度地图在海外地区，返回的是WGS84坐标。

谷歌地图和高德地图采用的是GCJ02地理坐标系。

验证了一下，谷歌地图和高德地图使用的经纬的确不一样，以天安门为例，高德地图上天安门的经纬度是116.397451,39.909187，而百度地图上该处经纬度是国家大剧院，二者离得很远的。

不同坐标系之间的相互转换

百度的地图API中支持坐标转换，提供将非百度坐标转换成百度坐标，但是根据相关法律规定，暂不支持将任何一种坐标系转换为WGS84类型(出自参考文献2)。

高德地图在高德开放平台->开发支持-> Web服务 API > 开发指南 > API文档 > 坐标转换中也提供了非高德坐标转换为高德坐标的方法。

据参考文献1所说，互联网地图在国内必须至少使用GCJ02进行首次加密，不允许直接使用WGS84坐标下的地理数据，同时任何坐标系均不可转换为WGS84坐标。

但是参考文献3中提供了三种坐标系相互转换的代码，试了试GCJ02转BD09，还是以天安门做例，高德的天安门坐标转过去正好是百度的天安门坐标，其他的暂时还没有试，下面是代码摘录：

import json
import urllib
import math

x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626  # π
a = 6378245.0  # 长半轴
ee = 0.00669342162296594323  # 偏心率平方

def gcj02_to_bd09(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]


def bd09_to_gcj02(bd_lon, bd_lat):
    """
    百度坐标系(BD-09)转火星坐标系(GCJ-02)
    百度——>谷歌、高德
    :param bd_lat:百度坐标纬度
    :param bd_lon:百度坐标经度
    :return:转换后的坐标列表形式
    """
    x = bd_lon - 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)
    gg_lng = z * math.cos(theta)
    gg_lat = z * math.sin(theta)
    return [gg_lng, gg_lat]


def wgs84_to_gcj02(lng, lat):
    """
    WGS84转GCJ02(火星坐标系)
    :param lng:WGS84坐标系的经度
    :param lat:WGS84坐标系的纬度
    :return:
    """
    if out_of_china(lng, lat):  # 判断是否在国内
        return [lng, lat]
    dlat = _transformlat(lng - 105.0, lat - 35.0)
    dlng = _transformlng(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]


def gcj02_to_wgs84(lng, lat):
    """
    GCJ02(火星坐标系)转GPS84
    :param lng:火星坐标系的经度
    :param lat:火星坐标系纬度
    :return:
    """
    if out_of_china(lng, lat):
        return [lng, lat]
    dlat = _transformlat(lng - 105.0, lat - 35.0)
    dlng = _transformlng(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]


def bd09_to_wgs84(bd_lon, bd_lat):
    lon, lat = bd09_to_gcj02(bd_lon, bd_lat)
    return gcj02_to_wgs84(lon, lat)


def wgs84_to_bd09(lon, lat):
    lon, lat = wgs84_to_gcj02(lon, lat)
    return gcj02_to_bd09(lon, lat)


def _transformlat(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


def _transformlng(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


def out_of_china(lng, lat):
    """
    判断是否在国内，不在国内不做偏移
    :param lng:
    :param lat:
    :return:
    """
    return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)

if __name__ == '__main__':
	print(gcj02_to_bd09(116.397451,39.909187))  # 把输出的坐标在百度地图上看了看，还是挺准的

计算两个经纬度之间的距离

python代码，取自参考文献4

EARTH_REDIUS = 6378.137

def rad(d):
    return d * pi / 180.0

def getDistance(lat1, lng1, lat2, lng2):
    radLat1 = rad(lat1)
    radLat2 = rad(lat2)
    a = radLat1 - radLat2
    b = rad(lng1) - rad(lng2)
    s = 2 * math.asin(math.sqrt(math.pow(sin(a/2), 2) + cos(radLat1) * cos(radLat2) * math.pow(sin(b/2), 2)))
    s = s * EARTH_REDIUS
    return s

返回单位是千米。

在高德地图上测过两个点的距离，跟这个计算得来的基本相符。

关于计算公式的意义和相关解析，可以参照参考文献5，虽然参考文献5是C#代码，不过它详细解释了计算的原理。

Java代码，参照参考文献6

public class DistanceUtil {
private static final double EARTH_RADIUS = 6378137;//赤道半径
private static double rad(double d){
return d * Math.PI / 180.0;
}

/ **
 * 经纬度获取距离
 * @param lon1 第一个经度
 * @param lat1 第一个纬度
 * @param lon2  
 * @param lat2
 * @return
 */

public static double getDistance(double lon1,double lat1,double lon2, double lat2) {
    double radLat1 = rad(lat1);
    double radLat2 = rad(lat2);
    double a = radLat1 - radLat2;
    double b = rad(lon1) - rad(lon2);
    double s = 2 *Math.asin(Math.sqrt(Math.pow(Math.sin(a/2),2)+Math.cos(radLat1)*Math.cos(radLat2)*Math.pow(Math.sin(b/2),2))); 
    s = s * EARTH_RADIUS;    
   return s/1000;//单位千米
}

public static void main(String[] args){  
    double distance = getDistance(自己的经度,自己的纬度,对方的经度,对方的纬度);  
    System.out.println("Distance is:"+distance);  
}  
}

各种辅助网站

国家统计局有提供，2018年统计用区划代码和城乡划分代码，这里提供到每个省市区镇(街道)村(居委会)的区划代码，通过这个也能得到详细的行政区划层级关系。比如说，通过这个我可以知道每个城市都有哪些村。除了这个官方网站之外，还有，行政区划这个网站也能达到类似效果。中华人民共和国民政部-行政区划代码 这里有每年的最新的行政区划代码，但是仅限于县以上。

参考文献

1. 坐标系说明书
2. 服务常见问题
3. 三种常见经纬度坐标系的转化
4. 计算两个经纬度之间的距离(python算法)
5. 通过经纬度坐标计算距离的方法（经纬度距离计算）
6. java根据经纬度获取距离