openlayers6【八】高德地图瓦片位置纠偏,将高德底图瓦片位置转换为EPSG:4326和EPSG:3857

369 阅读4分钟

实现功能

  1. 实现高德底图瓦片整体纠偏
  2. 使用两个图层来演示瓦片整体纠偏后的位置
import 'ol/ol.css'
import * as control from 'ol/control'
import { Feature, Map, View } from 'ol/index'
import { Tile as TileLayer, Vector as VectorLayer } from 'ol/layer'
import { Vector as VectorSource, XYZ } from 'ol/source'
import { Projection, addCoordinateTransforms } from 'ol/proj'
import { Fill, Stroke, Style, Circle as CircleStyle, Text } from 'ol/style'
import { Point } from 'ol/geom'

const forEachPoint = function(func) {
  return function(input, opt_output, opt_dimension) {
    const len = input.length
    const dimension = opt_dimension || 2
    let 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
  }
}
const gcj02 = {}
const PI = Math.PI
const AXIS = 6378245.0
const OFFSET = 0.00669342162296594323 // (a^2 - b^2) / a^2

function delta(wgLon, wgLat) {
  let dLat = transformLat(wgLon - 105.0, wgLat - 35.0)
  let dLon = transformLon(wgLon - 105.0, wgLat - 35.0)
  const radLat = wgLat / 180.0 * PI
  let magic = Math.sin(radLat)
  magic = 1 - OFFSET * magic * magic
  const 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) {
  let 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) {
  let 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 = {}

const RADIUS = 6378137
const MAX_LATITUDE = 85.0511287798
const RAD_PER_DEG = Math.PI / 180

sphericalMercator.forward = forEachPoint(function(input, output, offset) {
  const lat = Math.max(Math.min(MAX_LATITUDE, input[offset + 1]), -MAX_LATITUDE)
  const 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) {
  const output = gcj02.fromWGS84(input, opt_output, opt_dimension)
  return projzh.ll2smerc(output, output, opt_dimension)
}
projzh.gmerc2ll = function(input, opt_output, opt_dimension) {
  const 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

// 定义高德gcj02坐标系
const gcj02Extent = [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244]
const gcjMecator = new Projection({
  code: 'GCJ-02',
  extent: gcj02Extent,
  units: 'm'
})
// 注册GCJ02坐标系转换EPSG:4326
addCoordinateTransforms('EPSG:4326', gcjMecator, projzh.ll2gmerc, projzh.gmerc2ll)
// 注册GCJ02坐标系转换EPSG:3857
addCoordinateTransforms('EPSG:3857', gcjMecator, projzh.smerc2gmerc, projzh.gmerc2smerc)

var map
function initMap() {
  map = new Map({
    target: 'eguid-map',
    controls: control.defaults({ attribution: false, zoom: true, rotate: true }),
    view: new View({
      center: [121.499731, 31.2397],
      maxZoom: 18,
      zoom: 9,
      minZoom: 1,
      projection: 'EPSG:4326' // 可以支持EPSG:3857和EPSG:4326
    })
  })
  // 添加高德地图图层
  const xyzLayer = new TileLayer({
    source: new XYZ({
      url: 'http://wprd0{1-4}.is.autonavi.com/appmaptile?lang=zh_cn&size=1&style=7&x={x}&y={y}&z={z}',
      projection: gcjMecator // GCJ02扩展坐标系
    })
  })
  map.addLayer(xyzLayer)

  // 未纠偏前位置
  const vectorLayer = new VectorLayer({
    source: new VectorSource({
      features: [new Feature({
        type: 'point',
        geometry: new Point([121.499656, 31.239925])
      })]
    }),
    style: new Style({
      image: new CircleStyle({
        radius: 5,
        fill: new Fill({ color: 'black' }),
        stroke: new Stroke({
          color: 'white',
          width: 2
        })
      }),
      text: new Text({
        text: '底图纠偏后的高德原坐标',
        font: '14px Calibri,sans-serif',
        fill: new Fill({
          color: 'rgba(0, 0, 0, 1)'
        }),
        backgroundFill: new Fill({
          color: 'rgba(250, 250, 250, 0.9)'
        }),
        padding: [2, 2, 2, 2],
        textBaseline: 'bottom',
        offsetY: -12
      })
    })
  })

  map.addLayer(vectorLayer)

  // 使用GCJ02转换成wgs84坐标系后位置
  const vectorLayer2 = new VectorLayer({
    source: new VectorSource({
      features: [new Feature({
        type: 'point',
        geometry: new Point(gcj02towgs84(121.499656, 31.239925))
      })]
    }),
    projection: gcjMecator, // GCJ02扩展坐标系
    style: new Style({
      image: new CircleStyle({
        radius: 5,
        fill: new Fill({ color: 'red' }),
        stroke: new Stroke({
          color: 'white',
          width: 2
        })
      }),
      text: new Text({
        text: '底图纠偏后,高德坐标转换为wgs84坐标',
        font: '14px Calibri,sans-serif',
        fill: new Fill({
          color: 'rgba(0, 0, 0, 1)'
        }),
        backgroundFill: new Fill({
          color: 'rgba(250, 250, 250, 0.9)'
        }),
        padding: [2, 2, 2, 2],
        textBaseline: 'bottom',
        offsetY: -12
      })
    })
  })

  map.addLayer(vectorLayer2)
}

const a = 6378245.0
const ee = 0.00669342162296594323
/**
 * 经度转换
 * @param { Number } lng
 * @param { Number } lat
 */
function 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 * 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
}

/**
 * 纬度转换
 * @param { Number } lng
 * @param { Number } lat
 */
function 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 * 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
}

/**
 * GCJ02(高德)转WGS84(google)
 * @param { Number } lng:需要转换的经纬
 * @param { Number } lat:需要转换的纬度
 * @return { Array } result: 转换后的经纬度数组
 */
function gcj02towgs84(lng, lat) {
  if (outOfChina(lng, lat)) {
    return [lng, lat]
  } else {
    let dlat = transformlat(lng - 105.0, lat - 35.0)
    let dlng = transformlng(lng - 105.0, lat - 35.0)
    const radlat = lat / 180.0 * PI
    let magic = Math.sin(radlat)
    magic = 1 - ee * magic * magic
    const sqrtmagic = Math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * PI)
    dlng = (dlng * 180.0) / (a / sqrtmagic * Math.cos(radlat) * PI)
    const mglat = lat + dlat
    const mglng = lng + dlng
    return [lng * 2 - mglng, lat * 2 - mglat]
  }
}

测试效果演示

image.png

如此,我们便实现了自定义GCJ02坐标系,来实现对高德地图栅格/瓦片的纠偏!