MapBox 源码之墨卡托投影

MercatorCoordinate（墨卡托坐标）

MapBox 源码中 MercatorCoordinate (TileCoordinate) 是一种3D投影坐标。使用的是web墨卡托投影(EPSG:3857)，其单位略有不同:

• 通常墨卡托坐标范围是[-20037508，-20037508， 20037508， 20037508]，但是在Mapbox投影世界的最大宽度不是以米为单位，是按照1个作为整个投影世界的宽度， 范围为[0, 0, 1, 1]。
• 坐标轴的原点在西北角而不是中间， MercatorCoordinate(0,0,0)是mercator世界的西北角，而MercatorCoordinate(1,1,0)是东南角。

从源码中一下可以分析得出：

屏幕坐标转换为 MercatorCoordinate（墨卡托坐标）

// 屏幕坐标转换为经纬度坐标
LatLng TransformState::screenCoordinateToLatLng(const ScreenCoordinate& point, LatLng::WrapMode wrapMode) const {
     // 屏幕坐标转换为墨卡托投影坐标
     auto coord = screenCoordinateToTileCoordinate(point, 0);
     return Projection::unproject(coord.p, 1 / util::tileSize, wrapMode);
}

TileCoordinate TransformState::screenCoordinateToTileCoordinate(const ScreenCoordinate& point, uint8_t atZoom) const {
     #..............忽略其它无关代码..............
     mat4 mat = coordinatePointMatrix();

     mat4 inverted;
     bool err = matrix::invert(inverted, mat);

     double flippedY = size.height - point.y;

     vec4 coord0;
     vec4 coord1;
     vec4 point0 = {{ point.x, flippedY, 0, 1 }};
     vec4 point1 = {{ point.x, flippedY, 1, 1 }};

     //矩阵变换得到墨卡托投影坐标
     matrix::transformMat4(coord0, point0, inverted);
     matrix::transformMat4(coord1, point1, inverted);

     #..............忽略其它无关代码..............
     
     // / scale 严格按照 瓦片大小（512）来进行[0, 1] 换算
     Point<double> p = util::interpolate(p0, p1, t) / scale * static_cast<double>(1 << atZoom);   
     return { { p.x, p.y }, static_cast<double>(atZoom) };
}


mat4 TransformState::coordinatePointMatrix() const {
     mat4 proj;
     getProjMatrix(proj);
     matrix::scale(proj, proj, util::tileSize, util::tileSize, 1);
     matrix::multiply(proj, getPixelMatrix(), proj);
     return proj;
}

// NDC坐标归一化 
// 屏幕坐标系左上角为[0, 0] 右下递增  NDC坐标系以屏幕中心[0, 0], 右上递增
// 屏幕坐标系 * invert（getPixelMatrix）= NDC坐标系
mat4 TransformState::getPixelMatrix() const {
     mat4 m;
     matrix::identity(m);
     matrix::scale(m, m,
                   static_cast<double>(size.width) / 2, -static_cast<double>(size.height) / 2, 1);
     matrix::translate(m, m, 1, -1, 0);
     return m;
}

void TransformState::getProjMatrix(mat4& projMatrix, uint16_t nearZ, bool aligned) const {
     if (size.isEmpty()) {
         return;
     }
     #..............忽略其它无关代码..............
     const double dx = pixel_x() - size.width / 2.0f, dy = pixel_y() - size.height / 2.0f;
     
     // 中心点坐标偏移了瓦片坐标的一半, 移动到了左上角
     matrix::translate(projMatrix, projMatrix, dx, dy, 0);
     
     #..............忽略其它无关代码..............
     }
 }
 
double TransformState::pixel_x() const {
     const double center = (size.width - Projection::worldSize(scale)) / 2;
     return center + x;
}

double TransformState::pixel_y() const {
     const double center = (size.height - Projection::worldSize(scale)) / 2;
     return center + y;
}

static double worldSize(double scale) {
     return scale * util::tileSize;
}

MercatorCoordinate（墨卡托坐标）转换为经纬度坐标

MapBox采用的Web墨卡托投影，又叫球面墨卡托投影。接收的输入是 WGS84 的经纬度，但在投影时不再把地球当做椭球而当做半径为6378137米的标准球体，以简化计算。简化过的 Web墨卡托投影，球面坐标系经纬度与平面坐标系坐标之间的转换公式如下：

从以下MapBox中的源码分析可以得到：

constexpr double tileSize = 512;
constexpr int32_t EXTENT = 8192;

constexpr double DEG2RAD = M_PI / 180.0;
constexpr double RAD2DEG = 180.0 / M_PI;
constexpr double M2PI = M_PI * 2;
constexpr double EARTH_RADIUS_M = 6378137;
constexpr double LATITUDE_MAX = 85.051128779806604;
constexpr double LONGITUDE_MAX = 180;
constexpr double DEGREES_MAX = 360;


static LatLng unproject(const Point<double>& p, double scale, LatLng::WrapMode wrapMode = LatLng::Unwrapped) {
    auto p2 = p * util::DEGREES_MAX / worldSize(scale);
    return LatLng {
        util::DEGREES_MAX / M_PI * std::atan(std::exp((util::LONGITUDE_MAX - p2.y) * util::DEG2RAD)) - 90.0,
        p2.x - util::LONGITUDE_MAX,
        wrapMode
    };
}

static Point<double> project_(const LatLng& latLng, double worldSize) {
    const double latitude = util::clamp(latLng.latitude(), -util::LATITUDE_MAX, util::LATITUDE_MAX);
    return Point<double> {
        util::LONGITUDE_MAX + latLng.longitude(),
        util::LONGITUDE_MAX - util::RAD2DEG * std::log(std::tan(M_PI / 4 + latitude * M_PI / util::DEGREES_MAX))
    } * (worldSize / util::DEGREES_MAX);
}

经度转换

经度区间[-180, 180] 对应 [0, 1]，可简化为JS版本：

function mercatorXfromLng(lng: number): number { 
     return (180 + lng) / 360; 
} 
     
function lngFromMercatorX(x: number): number { 
     return x * 360 - 180; 
}

维度

维度区间[-180, 180] 对应 [0, 1]，可简化为JS版本：

function mercatorYfromLat(lat: number): number { 
    return (180 - (180 / Math.PI * Math.log(Math.tan(Math.PI / 4 + lat * Math.PI / 360)))) / 360; 
} 

function latFromMercatorY(y: number): number { 
    const y2 = 180 - y ; 
    return 360 / Math.PI * Math.atan(Math.exp(y2 * Math.PI / 180)) - 90; 
}

Refrence

• MapBoxGl Mercator Coordinate

• WEB地图基本原理：地图投影和坐标转换