关于openlayers引入高德地图出现坐标偏移的纠偏方法

之前通过openlayers成功将百度，高德地图引入系统，但是由于百度是BD-09坐标，高德是GCJ-02坐标系，所以当用WGS-84的坐标建立点，线，面时，这些要素就无法显示在正确的位置上，当然，这时候我们可以通过网上提供的坐标转换公式，对传过来的坐标进行转换，从而适应不同的地图，但是这种方法比较麻烦，每次都要转换，而且当要切换两种地图时候，还要对地图要素进行重新创建，所以一劳永逸的方法是对地图瓦片进行纠偏，刚好最近看到一篇博客https://blog.csdn.net/u013594477/article/details/83988055，实现了对百度地图坐标偏移的方法，并且亲测成功，所以照葫芦画瓢，就利用此方法实现对高德地图纠偏。

第一步：用openlayers加载高德在线地图

很简单，百度一搜一大堆，核心代码如下：

      new ol.layer.Tile({
        source:new ol.source.XYZ({
          url: 'http://wprd0{1-4}.is.autonavi.com/appmaptile?lang=zh_cn&size=1&style=7&x={x}&y={y}&z={z}'
        })
      }),

效果如下：

第二步：重新对切片进行投影

我们知道，目前openlayers只支持主流的坐标系，而对于GCJ-02以及BD-09是没有的，因此，我们要对瓦片进行纠偏。

1.定义GCJ-02投影

  const gcj02Extent = [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244];
  const gcjMecator = new ol.proj.Projection({
    code: "GCJ-02",
    extent: gcj02Extent,
    units: "m"
  });
  ol.proj.addProjection(gcjMecator);

2.进行坐标转换，将4326转gcjMecator，主要用ol.proj.addCoordinateTransforms方法，核心代码如下；

  ol.proj.addCoordinateTransforms("EPSG:4326", gcjMecator, projzh.ll2gmerc, projzh.gmerc2ll);
  ol.proj.addCoordinateTransforms("EPSG:3857", gcjMecator, projzh.smerc2gmerc, projzh.gmerc2smerc);

通过定义可知道，这两端代码主要是为了在4326和3857坐标系和新建的gcjMecator坐标系之间建立一种坐标转换关系，当我们用传统的wgs84坐标时，系统自动映射到gcjMecator坐标系，核心代码如下

  var forEachPoint = function(func) {
    return function(input, opt_output, opt_dimension) {
      var len = input.length;
      var dimension = opt_dimension ? opt_dimension : 2;
      var output;
      if (opt_output) {
        output = opt_output;
      } else {
        if (dimension !== 2) {
          output = input.slice();
        } else {
          output = new Array(len);
        }
      }
      for (var offset = 0; offset < len; offset += dimension) {
        func(input, output, offset)
      }
      return output;
    };
  };
  var gcj02 = {};
  var i=0;
  var PI = Math.PI;
  var AXIS = 6378245.0;
  var OFFSET = 0.00669342162296594323;  // (a^2 - b^2) / a^2

  function delta(wgLon, wgLat) {
    var dLat = transformLat(wgLon - 105.0, wgLat - 35.0);
    var dLon = transformLon(wgLon - 105.0, wgLat - 35.0);
    var radLat = wgLat / 180.0 * PI;
    var magic = Math.sin(radLat);
    magic = 1 - OFFSET * magic * magic;
    var sqrtMagic = Math.sqrt(magic);
    dLat = (dLat * 180.0) / ((AXIS * (1 - OFFSET)) / (magic * sqrtMagic) * PI);
    dLon = (dLon * 180.0) / (AXIS / sqrtMagic * Math.cos(radLat) * PI);
    return [dLon, dLat];
  }

  function outOfChina(lon, lat) {
    if (lon < 72.004 || lon > 137.8347) {
      return true;
    }
    if (lat < 0.8293 || lat > 55.8271) {
      return true;
    }
    return false;
  }

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

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

  gcj02.toWGS84 = forEachPoint(
    function(input, output, offset) {
      var lng = input[offset];
      var lat = input[offset + 1];
      if (!outOfChina(lng, lat)) {
        var deltaD = delta(lng, lat);
        lng = lng - deltaD[0];
        lat = lat - deltaD[1];
      }
      output[offset] = lng;
      output[offset + 1] = lat;
    });

  gcj02.fromWGS84 = forEachPoint(function(input, output, offset) {
    var lng = input[offset];
    var lat = input[offset + 1];
    if (!outOfChina(lng, lat)) {
      var deltaD = delta(lng, lat);
      lng = lng + deltaD[0];
      lat = lat + deltaD[1];
    }
    output[offset] = lng;
    output[offset + 1] = lat;
  });

  var sphericalMercator = {}

  var RADIUS = 6378137;
  var MAX_LATITUDE = 85.0511287798;
  var RAD_PER_DEG = Math.PI / 180;

  sphericalMercator.forward = forEachPoint(function(input, output, offset) {
    var lat = Math.max(Math.min(MAX_LATITUDE, input[offset + 1]), -MAX_LATITUDE);
    var sin = Math.sin(lat * RAD_PER_DEG);

    output[offset] = RADIUS * input[offset] * RAD_PER_DEG;
    output[offset + 1] = RADIUS * Math.log((1 + sin) / (1 - sin)) / 2;
  });

  sphericalMercator.inverse = forEachPoint(function(input, output, offset) {
    output[offset] = input[offset] / RADIUS / RAD_PER_DEG;
    output[offset + 1] = (2 * Math.atan(Math.exp(input[offset + 1] / RADIUS)) - (Math.PI / 2)) / RAD_PER_DEG;
  });


  var projzh = {};
  projzh.ll2gmerc = function(input, opt_output, opt_dimension) {
    let output =gcj02.fromWGS84(input, opt_output, opt_dimension);
    return projzh.ll2smerc(output, output, opt_dimension);
  };
  projzh.gmerc2ll = function(input, opt_output, opt_dimension) {
    let output = projzh.smerc2ll(input, input, opt_dimension);
    return gcj02.toWGS84(output, opt_output, opt_dimension);
  };
  projzh.smerc2gmerc = function(input, opt_output, opt_dimension) {
    let output = projzh.smerc2ll(input, input, opt_dimension);
    output = gcj02.fromWGS84(output, output, opt_dimension);
    return projzh.ll2smerc(output, output, opt_dimension);
  };
  projzh.gmerc2smerc = function(input, opt_output, opt_dimension) {
    let output = projzh.smerc2ll(input, input, opt_dimension);
    output = gcj02.toWGS84(output, output, opt_dimension);
    return projzh.ll2smerc(output, output, opt_dimension);
  };

  projzh.ll2smerc = sphericalMercator.forward;
  projzh.smerc2ll = sphericalMercator.inverse;

3.在source里面加入新定义的gcjMecator投影

      new ol.layer.Tile({
        source:new ol.source.XYZ({
          projection:gcjMecator,
          url: 'http://wprd0{1-4}.is.autonavi.com/appmaptile?lang=zh_cn&size=1&style=7&x={x}&y={y}&z={z}'
        })
      }),

第三步：利用WGS84坐标建立要素

至此，地图纠偏已完成，我们来用wgs84坐标验证一下，由于openstreetmap坐标是基于wgs84的，那么我们从地图中任意选取一个点记录一下。

我们建个点要素

  //创建点
  let saoguan = new ol.Feature({
    geometry:new ol.geom.Point(ol.proj.transform([117.2700, 31.8632],'EPSG:4326','EPSG:3857'))
  });

  saoguan.setStyle(new ol.style.Style({
      image:new ol.style.Icon({
        src:'//a.amap.com/jsapi_demos/static/demo-center/icons/poi-marker-default.png'
      })
    })
  );
  let source = new ol.source.Vector({
    features:[saoguan]
  });
  let ptlayer = new ol.layer.Vector({
    source: source
  });

结果显示地图位置为，可以看到位置已经实现了纠偏效果

总结：当然，虽然能过实现坐标纠偏，但是我们可以看到还是有些误差的，而且貌似地图会拉伸，有点形变，这些问题还有待继续研究……