三维green公式_进阶-Green-从零开始写GIS JS API（六）

最新推荐文章于 2022-09-03 17:54:08 发布

weixin_39934869

最新推荐文章于 2022-09-03 17:54:08 发布

阅读量466

点赞数

文章标签：三维green公式

本文链接：https://blog.csdn.net/weixin_39934869/article/details/112379690

版权

前言

本篇主题——坐标&高德地图（AMap）、百度地图（BMap）、谷歌地图（GMap）。坐标这个主题对于GIS是绕不开的，网络上充斥着不少尝试去描述和阐明坐标系的原理、坐标系的转换等等相关文章，其中有些是好的，有些是误人子弟。一直想写一篇关于坐标系基础的文章，但转念一想：一是很多内容确实在其它文章已提到，再写只是重复（本文会参考一些有关坐标系的文章，并给出链接）；二是坐标系其实很简单，简单到一句话就可以说清：就是矩阵变换，二维矩阵对应（平移和缩放）、三维矩阵对应（平移、缩放和旋转），没了，仅此；三是，很多概念只有动笔后才能出真知。

因此，希望本文，在实现高德地图（AMap）、百度地图（BMap）集成的过程中，再从另一角度来描述坐标系。

准备知识

在通读本文前，请先通过维基百科了解墨卡托（Mercartor）（EPSG：3359）以及Web墨卡托（Web-Mercator）（EPSG：3857）。

同时，可以了解下文：

国内主要地图瓦片坐标系定义及计算原理cntchen.github.io

上文重复内容，将不在本文赘述。

注：此文不适合初学者，初学者无法从本文中学习坐标系相关知识，请先了解坐标投影相关知识以及推荐的上文。

高德地图&谷歌地图

谷歌地图作为Web地图的鼻祖，几乎影响了所有其它Web地图的方方面面，从底层API的设计与实现，到表现层功能交互，可以这么说大多数网络地图从一开始都是参考谷歌地图。这说明谷歌地图设计的好。没错，但是说句不尊敬的话，谷歌地图API的设计与GIS一些概念真的不符，这一点远没有ArcGIS开发架构师理解透彻。好，扯远了。而这一开始其中又属高德地图跟的最紧，可以翻翻老版本的API。所以高德地图与谷歌地图，从投影，到切片方式，几乎一样，唯独的区别就是GCJ02的坐标偏移（谷歌中国地图，也有偏移，所以是完全一样。哦，配色和样式不一样。。。）。

讲坐标，离开参考离开图，几乎是天书，下文就列出一个全文的参考点，请记住：

五道口，华清嘉园GPS坐标：（116.327158，39.990912）

注意：该参考点在百度地图API用于说明坐标偏移，此处沿用。另要注意这是GPS坐标，即是GPS测量得到，非经过任何偏移。 再注意：本文所有Demo基于green-gis 0.0.5版本，该版本支持gcj02以及bd09。

首先，在正常Web-Mercator下，该坐标的显示：

红色箭头才是正确位置，为什么？这就是GCJ02（火星坐标系）要达到的目的。不管真实目的是什么，但明眼人一看就知道这几乎是掩耳盗铃。（WGS84 与 GCJ02 转换公式在Github或百度上一搜一大把。。。）

该例源码在：https://stackblitz.com/edit/green-gis-amap-gps

我们看下WebMercator的源码实现：（此处完全借用Leaflet的实现，借花献佛）

export class WebMercator extends Projection{
    static R: number = 6378137;
    //投影后的平面坐标范围
    get bound(): Bound {
        return new Bound(- Math.PI * WebMercator.R, Math.PI * WebMercator.R, Math.PI * WebMercator.R, -Math.PI * WebMercator.R);
    }
    //经纬度转平面坐标
    project([lng, lat]): number[] {
        //from leaflet & wiki
        const d = Math.PI / 180, sin = Math.sin(lat * d);
        return [WebMercator.R * lng * d,  WebMercator.R * Math.log((1 + sin) / (1 - sin)) / 2];
    }
    //平面坐标转经纬度
    unproject([x, y]): number[] {
        const d = 180 / Math.PI;
        return  [x * d / WebMercator.R, (2 * Math.atan(Math.exp(y / WebMercator.R)) - (Math.PI / 2)) * d];
    }
}

Leaflet的实现是最标准的，所以此处借用。需要注意的是Bound边界，具体为什么，请参考准备知识中的一文。此处做一伏笔，着重提醒的原因是，百度地图则完全不同！！

如何解决，参见GCJ02（此处转换参考https://github.com/wandergis/coordtransform）

export class GCJ02 extends Projection{
    static R = 6378137.0;
    static ee = 0.00669342162296594323;
    private _type: LatLngType;
    constructor(type: LatLngType = LatLngType.GPS) {
        super();
        this._type = type;
    }
    //投影后的平面坐标范围
    get bound(): Bound {
        return new Bound(- Math.PI * GCJ02.R, Math.PI * GCJ02.R, Math.PI * GCJ02.R, -Math.PI * GCJ02.R);
    }
    //经纬度转平面坐标
    project([lng, lat]): number[] {
        if (this._type == LatLngType.GPS) {
            [lng, lat] = GCJ02.wgs84togcj02(lng, lat);
        }
        //from leaflet & wiki
        const d = Math.PI / 180, sin = Math.sin(lat * d);
        return [GCJ02.R * lng * d,  GCJ02.R * Math.log((1 + sin) / (1 - sin)) / 2];
    }
    //平面坐标转经纬度
    unproject([x, y]): number[] {
        const d = 180 / Math.PI;
        return  [x * d / GCJ02.R, (2 * Math.atan(Math.exp(y / GCJ02.R)) - (Math.PI / 2)) * d];
    }

    //from https://github.com/wandergis/coordtransform
    /**
     * WGS-84 转 GCJ-02
     * @param lng
     * @param lat
     * @returns {*[]}
     */
    static wgs84togcj02(lng, lat) {
        var dlat = this._transformlat(lng - 105.0, lat - 35.0);
        var dlng = this._transformlng(lng - 105.0, lat - 35.0);
        var radlat = lat / 180.0 * Math.PI;
        var magic = Math.sin(radlat);
        magic = 1 - GCJ02.ee * magic * magic;
        var sqrtmagic = Math.sqrt(magic);
        dlat = (dlat * 180.0) / ((GCJ02.R * (1 - GCJ02.ee)) / (magic * sqrtmagic) * Math.PI);
        dlng = (dlng * 180.0) / (GCJ02.R / sqrtmagic * Math.cos(radlat) * Math.PI);
        var mglat = lat + dlat;
        var mglng = lng + dlng;
        return [mglng, mglat]
    };

    /**
     * GCJ-02 转换为 WGS-84
     * @param lng
     * @param lat
     * @returns {*[]}
     */
    static gcj02towgs84(lng, lat) {
        var dlat = this._transformlat(lng - 105.0, lat - 35.0);
        var dlng = this._transformlng(lng - 105.0, lat - 35.0);
        var radlat = lat / 180.0 * Math.PI;
        var magic = Math.sin(radlat);
        magic = 1 - GCJ02.ee * magic * magic;
        var sqrtmagic = Math.sqrt(magic);
        dlat = (dlat * 180.0) / ((GCJ02.R * (1 - GCJ02.ee)) / (magic * sqrtmagic) * Math.PI);
        dlng = (dlng * 180.0) / (GCJ02.R / sqrtmagic * Math.cos(radlat) * Math.PI);
        var mglat = lat + dlat;
        var mglng = lng + dlng;
        return [lng * 2 - mglng, lat * 2 - mglat]
    };

    static _transformlat(lng, lat) {
        var ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * Math.sqrt(Math.abs(lng));
        ret += (20.0 * Math.sin(6.0 * lng * Math.PI) + 20.0 * Math.sin(2.0 * lng * Math.PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(lat * Math.PI) + 40.0 * Math.sin(lat / 3.0 * Math.PI)) * 2.0 / 3.0;
        ret += (160.0 * Math.sin(lat / 12.0 * Math.PI) + 320 * Math.sin(lat * Math.PI / 30.0)) * 2.0 / 3.0;
        return ret
    };

    static _transformlng(lng, lat) {
        var ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * Math.sqrt(Math.abs(lng));
        ret += (20.0 * Math.sin(6.0 * lng * Math.PI) + 20.0 * Math.sin(2.0 * lng * Math.PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(lng * Math.PI) + 40.0 * Math.sin(lng / 3.0 * Math.PI)) * 2.0 / 3.0;
        ret += (150.0 * Math.sin(lng / 12.0 * Math.PI) + 300.0 * Math.sin(lng / 30.0 * Math.PI)) * 2.0 / 3.0;
        return ret
    };


    //此判断欠妥，暂不采用！
    /**
     * 判断是否在国内，不在国内则不做偏移
     * @param lng
     * @param lat
     * @returns {boolean}
     */
    static out_of_china(lng, lat) {
        // 纬度 3.86~53.55, 经度 73.66~135.05
        return !(lng > 73.66 && lng < 135.05 && lat > 3.86 && lat < 53.55);
    };
}

注意比较标准Web-Mercator 与 GCJ02，两个实现。唯一区别就是这一句：

        if (this._type == LatLngType.GPS) {
            [lng, lat] = GCJ02.wgs84togcj02(lng, lat);
        }

目的就是在经纬度转平面坐标之前，先进行偏移。

那采用该实现后，修改上例，参考https://stackblitz.com/edit/green-gis-amap-gcj02

注意其中的改动map.setProjection(new GCJ02(LatLngType.GPS));因为Map默认采用标准Web-Mercator。

另，对于GCJ02构造，传入LatLngType经纬度坐标的类型，该枚举分为：

export enum LatLngType {
    GPS = 0,           //Default
    GCJ02 = 1,         //Just For China, AMap aka GaoDe
    BD09 = 2           //Just For China, BaiduMap
}

该参数存在的原因是：有些经纬度坐标本身已偏移，如从高德地图中导出的数据，比如前文一直在使用的某市区县边界图（chongqing.json）,该数据本身已完成偏移，因此在project时，直接进行Web-Mercator标准变换。再看这一句：

        if (this._type == LatLngType.GPS) {
            [lng, lat] = GCJ02.wgs84togcj02(lng, lat);
        }

这一过程，对于百度地图也是类似，可参见下文描述百度地图实现后的总结。

谷歌集成示例：https://stackblitz.com/edit/green-gis-google，key现在要信用卡，所以：

小结，对于高德地图与谷歌地图集成实现（谷歌集成示例：），如果坐标不在国内，可直接采用标准的Web-Mercator实现；如果坐标在国内，根据提供坐标（经纬度）是否经过偏移来决定，是否在标准变换前进行偏移。

另，此处着重说明，GCJ02虽继承自Projection，但请注意，GCJ02并非为新的投影方式，其只不过是在投影前，先进行经纬度的偏移。此处为了调用者更为清楚区别是否要进行偏移，尤其针对国外坐标的使用者，如增加Web-Mercator构造参数（增加LatLngType），虽可去掉GCJ02该类的实现，但会造成调用者的困扰：标准的Web-Mercator，何来GCJ02？

故，特此说明。

百度地图

百度地图，从产品的角度来说，真的在国内是非常优秀，个人平常开车导航或地图查询一般用百度地图，但开发却用高德地图。为啥，说不出。有了GCJ02，再来个BD09，二次偏移，为了不一样而不一样？总之，偏移这事是真的鸡肋，掩耳盗铃。投影，百花齐放就忍了，有因为参考椭球体不同，有了需求而等角或等积；而偏移，百花齐放真的忍不了，除了降低效率，增加集成难度，真的能提高数据完全和保密？掩耳盗铃。BD09二次偏移还能忍，对于好多地方所谓的地方坐标系，打着数据完全和保密的旗号，做了偏移，真实目的就是为了钱，非常直白，忍无可忍。对于数据收费，可以理解，但通过偏移进行所谓的数据保护，进行重复多次取费，这些行为，无疑是建设智慧城市的噩梦！

言归正传，百度地图虽然也采用Web-Mercator投影，但是由于它的切片方式完全不同与高德和谷歌，因此，上文中提到的Bound边界也完全不同，而这会影响Canvas设置当前的变换矩阵。

注：此处摸索了一下午，最后终于理解那篇参考文章中的一句话：百度平面坐标系的坐标原点与百度瓦片坐标原点相同，以瓦片等级18级为基准，规定18级时百度平面坐标的一个单位等于屏幕上的一个像素

这句话体现到Bound边界和BD09的实现中：

export class BD09 extends Projection{
    //百度平面坐标系的坐标原点与百度瓦片坐标原点相同，以瓦片等级18级为基准，规定18级时百度平面坐标的一个单位等于屏幕上的一个像素
    static TOTAL_PIXELS = 256 * Math.pow(2, 18);
    private _type: LatLngType;
    constructor(type: LatLngType = LatLngType.GPS) {
        super();
        this._type = type;
    }
    //投影后的平面坐标范围
    get bound(): Bound {
        return new Bound(- BD09.TOTAL_PIXELS/2, BD09.TOTAL_PIXELS/2, BD09.TOTAL_PIXELS/2, -BD09.TOTAL_PIXELS/2);
    }
    //经纬度转平面坐标
    project([lng, lat]): number[] {
        //from leaflet & wiki
        if (this._type == LatLngType.GPS) {
            [lng, lat] = GCJ02.wgs84togcj02(lng, lat);
            [lng, lat] = BD09.gcj02tobd09(lng, lat);
        } else if (this._type == LatLngType.GCJ02) {
            [lng, lat] = BD09.gcj02tobd09(lng, lat);
        }
        const projection =  new BMap.MercatorProjection();
        const pixel = projection.lngLatToPoint(new BMap.Point(lng, lat));
        return [pixel.x, pixel.y];
        /*const d = Math.PI / 180, sin = Math.sin(lat * d);
        return [WebMercator.R * lng * d,  WebMercator.R * Math.log((1 + sin) / (1 - sin)) / 2];*/
    }
    //平面坐标转经纬度
    unproject([x, y]): number[] {
        const projection =  new BMap.MercatorProjection();
        const point = projection.pointToLngLat(new BMap.Pixel(x, y));
        return [point.lng, point.lat];
        /*const d = 180 / Math.PI;
        return  [x * d / WebMercator.R, (2 * Math.atan(Math.exp(y / WebMercator.R)) - (Math.PI / 2)) * d];*/
    }

    //from https://github.com/wandergis/coordtransform
    /**
     * 百度坐标系 (BD-09) 与 火星坐标系 (GCJ-02) 的转换
     * 即 百度 转 谷歌、高德
     * @param bd_lng
     * @param bd_lat
     * @returns {*[]}
     */
    static bd09togcj02(bd_lng, bd_lat) {
        var x = bd_lng - 0.0065;
        var y = bd_lat - 0.006;
        var z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * Math.PI * 3000.0 / 180.0);
        var theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * Math.PI * 3000.0 / 180.0);
        var gg_lng = z * Math.cos(theta);
        var gg_lat = z * Math.sin(theta);
        return [gg_lng, gg_lat]
    };

    /**
     * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换
     * 即 谷歌、高德 转 百度
     * @param lng
     * @param lat
     * @returns {*[]}
     */
    static gcj02tobd09(lng, lat) {
        var z = Math.sqrt(lng * lng + lat * lat) + 0.00002 * Math.sin(lat * Math.PI * 3000.0 / 180.0);
        var theta = Math.atan2(lat, lng) + 0.000003 * Math.cos(lng * Math.PI * 3000.0 / 180.0);
        var bd_lng = z * Math.cos(theta) + 0.0065;
        var bd_lat = z * Math.sin(theta) + 0.006;
        return [bd_lng, bd_lat]
    };

}

所以，高德和谷歌的Bound，很容易理解，就是地球的周长，而百度很特别，真的特别，在18级下平面单位与屏幕像素1：1，即256 * Math.pow(2, 18)，记住，如不特别指出，1个切片 = 256 *256 （pixel）。

另外，对应高德，请特别留意：

        if (this._type == LatLngType.GPS) {
            [lng, lat] = GCJ02.wgs84togcj02(lng, lat);
            [lng, lat] = BD09.gcj02tobd09(lng, lat);
        } else if (this._type == LatLngType.GCJ02) {
            [lng, lat] = BD09.gcj02tobd09(lng, lat);
        }

如果坐标是真实GPS，则先GCJ02偏移，再BD09偏移；如坐标是GCJ02，则只进行BD09偏移；如坐标是BD09，则不做处理；这一步完成后，再进行标准Web-Mercator变换。

再看，百度地图集成示例：https://stackblitz.com/edit/green-gis-baidu（您的密钥）

总结

汇总一个简单的流程决策：

以上，截止到0.0.5版本，目前支持集成高德、百度、谷歌等，当然腾讯地图等只要投影和切片方式类似高德和谷歌，也可参考高德和谷歌的集成示例。

weixin_39934869

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
三维green公式_进阶-Green-从零开始写GIS JS API（六）

前言本篇主题——坐标&高德地图（AMap）、百度地图（BMap）、谷歌地图（GMap）。坐标这个主题对于GIS是绕不开的，网络上充斥着不少尝试去描述和阐明坐标系的原理、坐标系的转换等等相关文章，其中有些是好的，有些是误人子弟。一直想写一篇关于坐标系基础的文章，但转念一想：一是很多内容确实在其它文章已提到，再写只是重复（本文会参考一些有关坐标系的文章，并给出链接）；二是坐标系其实很简单，简单...
复制链接

扫一扫