以下內(nèi)容轉(zhuǎn)載自木的樹的文章《Web地圖呈現(xiàn)原理》

作者：木的樹

鏈接：https://www.cnblogs.com/dojo-lzz/p/9250637.html

來源：博客園

著作權(quán)歸作者所有。商業(yè)轉(zhuǎn)載請聯(lián)系作者獲得授權(quán)，非商業(yè)轉(zhuǎn)載請注明出處劫流。

本篇內(nèi)容為大家揭開地圖呈現(xiàn)原理，大家可通過騰訊位置服務官網(wǎng)了解地圖并體驗榜配！

地圖投影

對于接觸互聯(lián)網(wǎng)地圖的同學來說，最開始接觸的恐怕就是坐標轉(zhuǎn)換的過程了吕晌。由于地球是個近似橢球的形狀芥牌，有各種各樣的橢球模型來模擬地球，最著名的也就是GPS系統(tǒng)使用的WGS84橢球了聂使。但是這些橢球體的坐標使用的是經(jīng)緯度壁拉，單位是角度。目前我們的地圖大多是二維平面上展示柏靶，使用角度為基礎來計算多有不便弃理，所以有眾多數(shù)學家提出各種不同的轉(zhuǎn)換方式來將經(jīng)緯度表示的位置轉(zhuǎn)換成平面坐標，這個轉(zhuǎn)換過程地圖學上成為投影屎蜓。投影的方式多種多樣痘昌，對我們做互聯(lián)網(wǎng)地圖的來說，最重要的就是墨卡托投影的變體——Web墨卡托投影炬转。我們先來看一下墨卡托投影的轉(zhuǎn)換過程

image

（以赤道本初子午線為原點）

投影完畢后的結(jié)果就是：

image

先不要頭疼數(shù)學公式辆苔，已經(jīng)有很多類庫做好了代碼實現(xiàn)，比如leaflet：

L.Projection.Mercator = {
    R: 6378137,
    R_MINOR: 6356752.314245179,

    bounds: L.bounds([-20037508.34279, -15496570.73972], [20037508.34279, 18764656.23138]),

    project: function (latlng) {
        var d = Math.PI / 180,
            r = this.R,
            y = latlng.lat * d,
            tmp = this.R_MINOR / r,
            e = Math.sqrt(1 - tmp * tmp),
            con = e * Math.sin(y);

        var ts = Math.tan(Math.PI / 4 - y / 2) / Math.pow((1 - con) / (1 + con), e / 2);
        y = -r * Math.log(Math.max(ts, 1E-10));

        return new L.Point(latlng.lng * d * r, y);
    },

    unproject: function (point) {
        var d = 180 / Math.PI,
            r = this.R,
            tmp = this.R_MINOR / r,
            e = Math.sqrt(1 - tmp * tmp),
            ts = Math.exp(-point.y / r),
            phi = Math.PI / 2 - 2 * Math.atan(ts);

        for (var i = 0, dphi = 0.1, con; i < 15 && Math.abs(dphi) > 1e-7; i++) {
            con = e * Math.sin(phi);
            con = Math.pow((1 - con) / (1 + con), e / 2);
            dphi = Math.PI / 2 - 2 * Math.atan(ts * con) - phi;
            phi += dphi;
        }

        return new L.LatLng(phi * d, point.x * d / r);
    }
};

接下來我們說一下互聯(lián)網(wǎng)地圖真正使用的投影——Web墨卡托或者也叫球形墨卡托扼劈。一般來說按照傳統(tǒng)地圖學的要求驻啤，一個投影坐標系都要有一個對應的橢球體，比如從WGS84的坐標轉(zhuǎn)換成國內(nèi)騰訊地圖或者百度地圖的坐標荐吵，都是要經(jīng)過一步橢球體轉(zhuǎn)換成gcj02橢球下的經(jīng)緯度然后才能打點骑冗。所以有沒有小伙伴在開發(fā)中使用Geolocation接口獲取的經(jīng)緯度直接傳入上面地圖api中打點發(fā)現(xiàn)誤差很大？就是因為沒有轉(zhuǎn)成gcj02橢球下的經(jīng)緯度先煎。但是web墨卡托這個投影其實并不符合地圖學的要求贼涩，它沒有對應的橢球體，它是谷歌自己造出來的（因為簡單）薯蝎，也可以說對任何橢球體都適用遥倦，但這種時候我們在表達一個位置信息時嚴格來說應當這樣表達：橢球下的Web墨卡托投影坐標是。

好了現(xiàn)在來說一下web墨卡托的轉(zhuǎn)換方式：

/*
 * @namespace Projection
 * @projection L.Projection.SphericalMercator
 *
 * Spherical Mercator projection — the most common projection for online maps,
 * used by almost all free and commercial tile providers. Assumes that Earth is
 * a sphere. Used by the `EPSG:3857` CRS.
 */

L.Projection.SphericalMercator = {

    R: 6378137,
    MAX_LATITUDE: 85.0511287798,

    project: function (latlng) {
        var d = Math.PI / 180,
            max = this.MAX_LATITUDE,
            lat = Math.max(Math.min(max, latlng.lat), -max),
            sin = Math.sin(lat * d);

        return new L.Point(
                this.R * latlng.lng * d,
                this.R * Math.log((1 + sin) / (1 - sin)) / 2);
    },

    unproject: function (point) {
        var d = 180 / Math.PI;

        return new L.LatLng(
            (2 * Math.atan(Math.exp(point.y / this.R)) - (Math.PI / 2)) * d,
            point.x * d / this.R);
    },

    bounds: (function () {
        var d = 6378137 * Math.PI;
        return L.bounds([-d, -d], [d, d]);
    })()
};

相對來說這個計算方法簡單一些占锯，但是這種投影有一些缺點袒哥，比如維度投影范圍只能在-85~85之間，一般來說沒什么關系烟央，反正一般人不會閑的沒事跑南北極去统诺。同時由于南北極特殊的位置，通常有一些專門的地圖來表達疑俭。

地圖的組織方式

可以觀察一下騰訊地圖在展示時粮呢，通常是一個正方形一個正方形的出現(xiàn)，這些正方形地圖上成為瓦片钞艇。下面我們來說一下地圖的組織方式啄寡。

如果地圖數(shù)據(jù)有一個G，那么在端上展示地圖時哩照，是把整個地圖數(shù)據(jù)全部下載下來好還是只把我們需要看的一部分下來好呢挺物。答案肯定是后者。那么有來一個問題飘弧，是每次都根據(jù)位置實時計算好還是提前將地圖分割成塊识藤，根據(jù)范圍加載瓦片好呢砚著？這個問題的答案不完全是瓦片，但絕大多數(shù)都是痴昧。所以現(xiàn)在互聯(lián)網(wǎng)地圖的組織形式就是金字塔結(jié)構(gòu)的瓦片地圖稽穆。

瓦片地圖金字塔模型是一種多分辨率層次模型，從瓦片金字塔的底層到頂層赶撰，分辨率越來越低舌镶，但表示的地理范圍不變（這張圖瓦片坐標是從左上角開始，通常谷歌系標準的瓦片是從左下角起始的）豪娜。

image

image

//示例來自：http://www.cnblogs.com/naaoveGIS/
//這里的像素是設備像素

現(xiàn)在假設地圖的坐標單位是米，dpi為96 惕虑；
    1英寸=2.54厘米坟冲； 
    1英寸=96像素； 
    最終換算的單位是米溃蔫； 
    如果當前地圖比例尺為1:125000000健提，則代表圖上1米等于實地125000000米； 
    米和像素間的換算公式： 
    1英寸=0.0254米=96像素 
    1像素=0.0254/96 米 
    則根據(jù)1：125000000比例尺伟叛，圖上1像素代表實地距離是 125000000*0.0254/96 = 33072.9166666667米私痹。

上面這個例子同樣可以由分辨率換算出比例尺。所以在互聯(lián)網(wǎng)地圖中我們先不去關心比例尺统刮，只關心分辨率的概念紊遵，假設瓦片的大小為256像素，每一級別的瓦塊數(shù)目為2^n * 2^n侥蒙。分辨率計算公式為：

r=6378137
resolution=2*PI*r/(2^zoom*256)

r為Web墨卡托投影時選取的地球半徑暗膜，2PIr代表地球周長，地球周長除以該級別下所有瓦片像素得到分辨率鞭衩。

注意這里的像素實際并不是設備像素学搜，而是一種參照像素（web中的css像素或者是安卓上的設備無關像素），比如在某些高清屏下（window.devicePixelRatio = 3）论衍，一參照像素寬度等于3設備像素瑞佩，這時候可能實際瓦片大小是512設備像素的，但是在顯示時仍然要把它顯示成256參照像素（css像素）坯台。

從經(jīng)緯度到地圖瓦片

現(xiàn)在要進入關鍵的一步炬丸！如何給定經(jīng)緯度來找出對應瓦片。這時候又要經(jīng)過幾個坐標轉(zhuǎn)換的過程：

1蜒蕾、經(jīng)緯度轉(zhuǎn)Web墨卡托稠炬；

2焕阿、Web墨卡托轉(zhuǎn)世界平面點；

3酸纲、世界平面點轉(zhuǎn)瓦片像素坐標捣鲸；

4、瓦片像素坐標轉(zhuǎn)瓦片行列號

我們再來看一下這些瓦片的組織方式闽坡，

image

可以看到這里的起始點是從左上角開始的，而經(jīng)緯度和Web墨卡托的起始點是赤道和本初子午線愁溜，在瓦片像素坐標的中心它的坐標是（256 * 2^n / 2, 256 * 2^n / 2）疾嗅，所以中間就有了世界平面點這一步，它是一個中間轉(zhuǎn)換的過程冕象。

所以上文中我們給出了經(jīng)緯度轉(zhuǎn)Web墨卡托的代碼代承。那么接下來就要把Web墨卡托坐標轉(zhuǎn)換成為世界平面點，這個坐標系的原點在左上角（0, 0）渐扮，在leaft中它認為這個坐標的原點在左上角為（0,0）论悴，坐標范圍為0~1。對應代碼：

// @method latLngToPoint(latlng: LatLng, zoom: Number): Point
    // Projects geographical coordinates into pixel coordinates for a given zoom.
    latLngToPoint: function (latlng, zoom) {
        var projectedPoint = this.projection.project(latlng),
            scale = this.scale(zoom);

        return this.transformation._transform(projectedPoint, scale);
    },
// @method scale(zoom: Number): Number
    // Returns the scale used when transforming projected coordinates into
    // pixel coordinates for a particular zoom. For example, it returns
    // `256 * 2^zoom` for Mercator-based CRS.
    scale: function (zoom) {
        return 256 * Math.pow(2, zoom);
    },

transform對應的代碼為：

/*
 * @class Transformation
 * @aka L.Transformation
 *
 * Represents an affine transformation: a set of coefficients `a`, `b`, `c`, `d`
 * for transforming a point of a form `(x, y)` into `(a*x + b, c*y + d)` and doing
 * the reverse. Used by Leaflet in its projections code.
 *
 * @example
 *
 * ```js
 * var transformation = new L.Transformation(2, 5, -1, 10),
 *     p = L.point(1, 2),
 *     p2 = transformation.transform(p), //  L.point(7, 8)
 *     p3 = transformation.untransform(p2); //  L.point(1, 2)
 * ```
 */


// factory new L.Transformation(a: Number, b: Number, c: Number, d: Number)
// Creates a `Transformation` object with the given coefficients.
L.Transformation = function (a, b, c, d) {
    this._a = a;
    this._b = b;
    this._c = c;
    this._d = d;
};

L.Transformation.prototype = {
    // @method transform(point: Point, scale?: Number): Point
    // Returns a transformed point, optionally multiplied by the given scale.
    // Only accepts actual `L.Point` instances, not arrays.
    transform: function (point, scale) { // (Point, Number) -> Point
        return this._transform(point.clone(), scale);
    },

    // destructive transform (faster)
    _transform: function (point, scale) {
        scale = scale || 1;
        point.x = scale * (this._a * point.x + this._b);
        point.y = scale * (this._c * point.y + this._d);
        return point;
    },

    // @method untransform(point: Point, scale?: Number): Point
    // Returns the reverse transformation of the given point, optionally divided
    // by the given scale. Only accepts actual `L.Point` instances, not arrays.
    untransform: function (point, scale) {
        scale = scale || 1;
        return new L.Point(
                (point.x / scale - this._b) / this._a,
                (point.y / scale - this._d) / this._c);
    }
};

對于Web墨卡托來說墓律，transformation的四個參數(shù)為：

/*
 * @namespace CRS
 * @crs L.CRS.EPSG3857
 *
 * The most common CRS for online maps, used by almost all free and commercial
 * tile providers. Uses Spherical Mercator projection. Set in by default in
 * Map's `crs` option.
 */

L.CRS.EPSG3857 = L.extend({}, L.CRS.Earth, {
    code: 'EPSG:3857',
    projection: L.Projection.SphericalMercator,

    transformation: (function () {
        var scale = 0.5 / (Math.PI * L.Projection.SphericalMercator.R);
        return new L.Transformation(scale, 0.5, -scale, 0.5);
    }())
});

L.CRS.EPSG900913 = L.extend({}, L.CRS.EPSG3857, {
    code: 'EPSG:900913'
});

這里要解釋一下在Web墨卡托中transform這四個參數(shù)的意思：scale代表球的周長分之一膀估，b和d都是0.5這代表赤道和本初子午線的交點在世界平面點的位置為(0.5, 0.5)；this._a * point.x + this._b代表x軸方向墨卡托坐標在世界平面點位置耻讽，c=-scale察纯，結(jié)合 this._c * point.y + this._d，能計算出y軸方向墨卡托在世界平面點位置针肥。至于c為什么是付的饼记，結(jié)合一下維度坐標的性質(zhì)，以上為正下為負慰枕，到了世界平面坐標中具则，負的緯度坐標要大于0.5。

接下來的兩步就比較簡單了具帮，世界平面點坐標轉(zhuǎn)像素坐標博肋，只要乘以各個軸上對應的像素長度：

256 * 2^zoom * coord.x
256 * 2^zoom * coord.y

在leaft中其實已經(jīng)在Transformation的_transform函數(shù)中坐了這一步。

剩下的我們已經(jīng)知道了像素坐標匕坯，就很容易求出對應的瓦片：

//tileSize = 256
xIndex = pixelCoord.x / tileSize;
yIndex = pixelCoord.y / tileSize;

注意谷歌系的瓦片都是以左下角為瓦片索引的起始束昵，所以對應的y方向上的瓦片計算方式為：

Math.pow(2, mapZoom) - yIndex - 1

加載一屏地圖

一般來說在實例化一個地圖時，都會給給Map構(gòu)造函數(shù)傳入一個zoom和一個center參數(shù)葛峻，3d情況下還會傳入theta和phi代表左右旋轉(zhuǎn)和上下翻轉(zhuǎn)锹雏，然后就會加載出一幅地圖。為了簡單起見我們先說說2D情況术奖，以leaflet為例

image

要加載一幅地圖礁遵，我們只需要知道屏幕四個點的經(jīng)緯度所在范圍內(nèi)的瓦片轻绞，再將這些瓦片按照一定的偏移坐標布置即可。

上面?zhèn)魅氲腸enter代表當前范圍的中心點佣耐，同時也是屏幕的中心點政勃，那么就可以求出該經(jīng)緯度對應的像素坐標，這個像素坐標就是屏幕中心點對應的瓦片像素坐標兼砖。這里的像素與我們的css像素一一對應奸远，利用屏幕范圍可得出屏幕四個角點的瓦片像素坐標。利用這四個點的瓦片坐標讽挟，可以求出當前屏幕的瓦片索引范圍懒叛，加載這些瓦片。

_getTiledPixelBounds: function (center) {
        var map = this._map,
            mapZoom = map._animatingZoom ? Math.max(map._animateToZoom, map.getZoom()) : map.getZoom(),
            scale = map.getZoomScale(mapZoom, this._tileZoom),
            pixelCenter = map.project(center, this._tileZoom).floor(),
            halfSize = map.getSize().divideBy(scale * 2);

        return new L.Bounds(pixelCenter.subtract(halfSize), pixelCenter.add(halfSize));
    },

_pxBoundsToTileRange: function (bounds) {
        var tileSize = this.getTileSize();
        return new L.Bounds(
            bounds.min.unscaleBy(tileSize).floor(),
            bounds.max.unscaleBy(tileSize).ceil().subtract([1, 1]));
    },

接下來要注意一些耽梅，這時候這些瓦片的坐標范圍肯定是大于屏幕的坐標范圍薛窥，所以要對所有的瓦片做一些偏移。計算過程比較簡單眼姐，屏幕坐標以左上點為原點诅迷，這個點對應的像素坐標已知，只要求出每個瓦片的左上角點的瓦片像素坐標與屏幕左上點的瓦片像素坐標做差值众旗，即可得出在css中的position的偏移值（高級點的用css3的translate）罢杉。下面我們可以看看leaflet的處理過程：

_setView: function (center, zoom, noPrune, noUpdate) {
        var tileZoom = Math.round(zoom);
        if ((this.options.maxZoom !== undefined && tileZoom > this.options.maxZoom) ||
            (this.options.minZoom !== undefined && tileZoom < this.options.minZoom)) {
            tileZoom = undefined;
        }

        var tileZoomChanged = this.options.updateWhenZooming && (tileZoom !== this._tileZoom);

        if (!noUpdate || tileZoomChanged) {

            this._tileZoom = tileZoom;

            if (this._abortLoading) { // 如果zoom要發(fā)生變化，停止當前所有tiles的加載逝钥；通過更改他們的onload onerror事件實現(xiàn)
                this._abortLoading();
            }

            // 1屑那、創(chuàng)建該級別容器瓦片
            // 2、 設置zIndex
            // 3艘款、設置本級別圖層基準點origin
            // 4持际、設置值本級別容器的偏移量
            this._updateLevels();
            // 1、得到世界的像素bounds
            // 2哗咆、得通過像素范圍除以tileSize得到能夠覆蓋世界的瓦片范圍
            // 3蜘欲、得到坐標系經(jīng)度和緯度范圍內(nèi)的瓦片范圍
            this._resetGrid();

            if (tileZoom !== undefined) {
                // 加載可視范圍內(nèi)的瓦片
                // 1、計算可視區(qū)域的像素范圍
                // 2晌柬、 將像素范圍轉(zhuǎn)變成瓦片格網(wǎng)范圍
                // 3姥份、計算一個buffer的格網(wǎng)范圍
                // 4、將不再當前范圍內(nèi)已加載的瓦片打上標簽
                // 5年碘、如果zoom發(fā)生變化重新進行setView
                // 6澈歉、將格網(wǎng)范圍內(nèi)的tile放入一個數(shù)組中
                // 7、對數(shù)組進行排序屿衅，靠近中心點的先加載
                // 8埃难、創(chuàng)建瓦片
                //     (1) 計算瓦片在地圖上的偏移量 coords * tileSize - origin
                //     (2) 加載瓦片數(shù)據(jù)（圖片或者矢量數(shù)據(jù)）
                //  (3) 設置圖片位置 setPosition
                this._update(center);
            }

            if (!noPrune) {
                this._pruneTiles(); // 移除不在范圍內(nèi)的tile; retainParent部分尚沒看懂，可能是按照瓦片金字塔保留
            }

            // Flag to prevent _updateOpacity from pruning tiles during
            // a zoom anim or a pinch gesture
            this._noPrune = !!noPrune;
        }
//將地圖的新中心點移到地圖中央
        this._setZoomTransforms(center, zoom);
    },

3D地圖的加載

互聯(lián)網(wǎng)地圖發(fā)展到現(xiàn)在出現(xiàn)了不少3D地圖，3D的計算過程有些復雜涡尘，尤其是設置了旋轉(zhuǎn)和俯視角度之后忍弛。不過我們可以從簡單情況說起。

3D地圖其實比2D多了一個環(huán)節(jié)考抄，那就是墨卡托與3D世界坐標细疚，3D世界與屏幕像素之間的轉(zhuǎn)換。如果我們不想自找麻煩川梅，通常3D坐標都是以米為單位疯兼，保持跟墨卡托的坐標單位一致，但是一般不直接以墨卡托的原點做3D世界的原點挑势，因為墨卡托坐標比較大镇防，后續(xù)計算精度是個問題。所以一般以用戶設置的center轉(zhuǎn)換成的墨卡托坐標為原點來建立3D的世界坐標系潮饱。

一般來講大家使用的都是透視投影，由于3D世界在屏幕上的投影時非線性的诫给，所以3D世界與屏幕像素之間的比值并不是固定的香拉。但一般對地圖來講，不考旋轉(zhuǎn)俯視情況下中狂，相機的視線軸與水平面是垂直關系凫碌，那么就可以利用相機投影面高度與屏幕高的比值求出3D世界單位與像素的比值，這個分辨率我們成為resolution2

_getPixelMeterRatio(target) {
    target = target ? target : this.controls.target;
    let distance = this.camera.position.distanceTo(target);
    let top = this.camera instanceof PerspectiveCamera ?
      (Math.tan(this.camera.fov / 2 * DEG2RAD) * distance) :
      this.camera.top / this.camera.zoom;
    let meterPerPixel = 2 * top / this.container.clientHeight;

    return meterPerPixel;
  }

前面章節(jié)中我們有一個地圖的瓦片像素分辨率resolution胃榕，只要讓這兩個分辨率的值相等盛险，就能計算出相機應當距離水平面的合適高度，將css像素與瓦片像素一比一對應起來勋又。但是這個時候不建議像2D那樣用四個點的瓦片像素坐標來計算瓦片索引苦掘，建議將屏幕四個點轉(zhuǎn)成3D坐標，3D坐標轉(zhuǎn)成墨卡托楔壤，墨卡托轉(zhuǎn)瓦片像素坐標鹤啡，然后再求瓦片索引。為什么要多此一舉蹲嚣，原因在于當俯視角度存在時递瑰，瓦片分辨率與resolution2并不相同，這時候的視野范圍是個梯形隙畜，但是我們可以將屏幕坐標轉(zhuǎn)成墨卡托坐標再來計算這個過程抖部。

image

image

是的沒錯，那個梯形就是你的手機屏幕议惰！至于這個梯形的計算過程慎颗，可以利用相機的fov、near、aspect等屬性計算出四條射線哗总，這四條射線與水平面的交點構(gòu)成了一個梯形几颜，這個梯形可以求出一個外包矩形，利用這個外包矩形求出瓦片的索引范圍讯屈。像mapbox中計算的方法相對巧妙一些蛋哭，沒有直接通過相機坐標系求射線方式，而是利用屏幕坐標求出ndc坐標涮母，通過將兩個ndc坐標的z值分別設置為0和1求出一條射線谆趾，然后將ndc坐標轉(zhuǎn)換成3d坐標，利用線性關系求出水平面的點（z=0）叛本。從而求出那個梯形沪蓬。

/**
     * Return all coordinates that could cover this transform for a covering
     * zoom level.
     * @param {Object} options
     * @param {number} options.tileSize
     * @param {number} options.minzoom
     * @param {number} options.maxzoom
     * @param {boolean} options.roundZoom
     * @param {boolean} options.reparseOverscaled
     * @param {boolean} options.renderWorldCopies
     * @returns {Array<Tile>} tiles
     */
    coveringTiles(
        options: {
            tileSize: number,
            minzoom?: number,
            maxzoom?: number,
            roundZoom?: boolean,
            reparseOverscaled?: boolean,
            renderWorldCopies?: boolean
        }
    ) {
        let z = this.coveringZoomLevel(options);
        const actualZ = z;

        if (options.minzoom !== undefined && z < options.minzoom) return [];
        if (options.maxzoom !== undefined && z > options.maxzoom) z = options.maxzoom;

        const centerCoord = this.pointCoordinate(this.centerPoint, z);
        const centerPoint = new Point(centerCoord.column - 0.5, centerCoord.row - 0.5);
// 利用屏幕四個點求ndc坐標，ndc坐標轉(zhuǎn)3D坐標来候，根據(jù)線性關系求出交點
        const cornerCoords = [
            this.pointCoordinate(new Point(0, 0), z),
            this.pointCoordinate(new Point(this.width, 0), z),
            this.pointCoordinate(new Point(this.width, this.height), z),
            this.pointCoordinate(new Point(0, this.height), z)
        ];
        return tileCover(z, cornerCoords, options.reparseOverscaled ? actualZ : z, this._renderWorldCopies)
            .sort((a, b) => centerPoint.dist(a.canonical) - centerPoint.dist(b.canonical));
    }

下面這個函數(shù)就是mapbox中求交點步驟的巧妙之處

pointCoordinate(p: Point, zoom?: number) {
        if (zoom === undefined) zoom = this.tileZoom;

        const targetZ = 0;
        // since we don't know the correct projected z value for the point,
        // unproject two points to get a line and then find the point on that
        // line with z=0

        const coord0 = [p.x, p.y, 0, 1];
        const coord1 = [p.x, p.y, 1, 1];

        vec4.transformMat4(coord0, coord0, this.pixelMatrixInverse);
        vec4.transformMat4(coord1, coord1, this.pixelMatrixInverse);

        const w0 = coord0[3];
        const w1 = coord1[3];
        const x0 = coord0[0] / w0;
        const x1 = coord1[0] / w1;
        const y0 = coord0[1] / w0;
        const y1 = coord1[1] / w1;
        const z0 = coord0[2] / w0;
        const z1 = coord1[2] / w1;

        const t = z0 === z1 ? 0 : (targetZ - z0) / (z1 - z0);

        return new Coordinate(
            interp(x0, x1, t) / this.tileSize,
            interp(y0, y1, t) / this.tileSize,
            this.zoom)._zoomTo(zoom);
    }