🛵 👊🏽 🎅🏼 从栅格地图创建图块 🎞️ 🦏 🛋️

不知何故，我对创建适合在OsmAnd和OpenLayers中使用的地图的问题感到困惑。那时我完全不了解GIS，因此我从头开始处理所有事情。

在本文中，我将向您介绍“研究”的结果，并提出一种算法，将任意栅格地图转换为应用程序可以理解的图块，并沿途熟悉诸如椭球，基准面，坐标系，投影等概念。

我们的地球不是球形的，甚至不是椭圆形的；这个复杂的数字称为大地水准面。事实是地球内部的质量分布不均匀，因此在某些地方地球略微凹陷，而在另一些地方则略微凸起。如果我们取一个单独的国家或大陆的领土，则如果该椭圆体的中心沿三个坐标轴相对于地球质心略有偏移，则其表面可以与某些椭圆体的表面对齐。这样的椭球称为参考椭球，它仅适用于描述为其创建地球的局部区域。与这个站点相距很远时，它与地球表面可能有很大的差异。中心与地球质心重合的椭球称为普通陆地椭球。明确，认为参考椭球比一般陆地更好地描述了地球的局部，但是一般陆地适用于地球的整个表面。

为了描述椭球，仅两个独立的值就足够了：赤道半径（通常用a表示）和极半径（b），但是通常使用极收缩率f =（ab）/ a代替第二个独立值。这是我们算法中需要作为对象的第一件事-椭球。对于不同年份的地球不同地区，不同的研究人员已经计算出许多参考椭球，有关它们的信息以数据的形式给出：a（以米为单位）和1 / f（无量纲）。奇怪的是，对于常见的陆地椭球体，也有许多不同的变体（不同的a，f），但差异不是很大，主要是由于确定a和f的方法不同。

struct Ellipsoid {
    char *name;
    double a;  /*  ()       */
    double b;  /*  ()               */
    double al; /*  (a-b)/a                        */
    double e2; /*   (a^2-b^2)/a^2 */
};

struct Ellipsoid Ellipsoid_WGS84 = {
    .name = "WGS84",
    .a  = 6378137.0,
    .al = 1.0 / 298.257223563,
};

struct Ellipsoid Ellipsoid_Krasovsky = {
    .name = "Krasovsky",
    .a  = 6378245.0,
    .al = 1.0 / 298.3,
};

该示例显示了两个椭球：用于卫星导航的普通地面WGS84和适用于欧洲和亚洲的Krasovsky参考椭球。

考虑另一个有趣的事实，事实是地球的形状很慢，但是却在变化，所以今天已经很好地描述了地表的椭球在一百年之内可能离现实还很远。这与普通的地球椭球无关，因为误差在同一误差范围内，但与参考椭球有关。在这里，我们得出另一个概念-基准。基准是椭球（a，f）的一组参数，其在地球内部的位移（对于参考椭球）在特定时刻固定。更准确地说，该基准可能不一定描述椭圆体，它可以是更复杂图形的参数，例如拟类曲面。

当然已经出现了一个问题：如何从一个椭球或基准移动到另一个？为此，每个椭球都必须具有一个大地坐标系：纬度和经度（phi，lambda），其转换是通过将坐标从一个坐标系转换到另一个坐标系来进行的。

有多种用于转换坐标的公式。您可以先将一个坐标系中的大地坐标转换为三维坐标X，Y，Z，并对其进行平移和转弯，然后将所得的三维坐标转换为另一坐标系中的大地坐标。您可以直接进行。因为由于所有公式都是无限收敛的级数，因此通常仅限于级数的几个成员以实现所需的精度。作为示例，我们将使用Helmert变换，这些变换是使用过渡到三维坐标的变换，它们由上述三个阶段组成。对于变换，除了两个椭圆形之外，还需要其他7个参数：沿三个轴的三个偏移，三个旋转角度和一个比例因子。您可能会猜到，所有参数都可以从基准中提取。但是在算法中，我们不会将这样的对象用作基准，而是将引入一个从一个坐标系到另一个坐标系的过渡对象，该过渡对象将包含：对两个椭圆体的引用和7个变换参数。这将是我们算法的第二个目标。

struct HelmertParam {
    char *src, *dst;
    struct Ellipsoid *esp;
    struct Ellipsoid *edp;
    struct {
        double dx, dy, dz;
        double wx, wy, wz;
        double ms;
    } p;
    //  
    double a,  da;
    double e2, de2;
    double de2__2, dxe2__2;
    double n, n__2e2;
    double wx_2e2__ro, wy_2e2__ro;
    double wx_n__ro, wy_n__ro;
    double wz__ro, ms_e2;
};

struct HelmertParam Helmert_SK42_WGS84 = {
    .src = "SK42",
    .dst = "WGS84",
    .esp = &Ellipsoid_Krasovsky,
    .edp = &Ellipsoid_WGS84,
    // SK42->PZ90->WGS84 (  51794-2001)
    .p = {23.92, -141.27, -80.9, 0, -0.35, -0.82, -0.12e-6},
};

该示例显示了用于从Pulkovo 1942坐标系转换为WGS84坐标系的参数。变换参数本身是一个单独的主题，对于某些坐标系来说是开放的，对于另一些则是凭经验选择的，因此它们的值在不同的来源中可能会略有不同。

除了参数外，还需要一个转换函数，它可以是直接的，而对于相反方向的转换，我们只需要相反方向的转换即可。我将跳过大量的Matan，并提供优化的功能。

void setupHelmert(struct HelmertParam *pp) {
    pp->a = (pp->edp->a + pp->esp->a) / 2;
    pp->da = pp->edp->a - pp->esp->a;
    pp->e2 = (pp->edp->e2 + pp->esp->e2) / 2;
    pp->de2 = pp->edp->e2 - pp->esp->e2;
    pp->de2__2 = pp->de2 / 2;
    pp->dxe2__2 = pp->de2__2 + pp->e2 * pp->da / pp->a;
    pp->n = 1 - pp->e2;
    pp->n__2e2 = pp->n / pp->e2 / 2;
    pp->wx_2e2__ro = pp->p.wx * pp->e2 * 2 * rad(0,0,1);
    pp->wy_2e2__ro = pp->p.wy * pp->e2 * 2 * rad(0,0,1);
    pp->wx_n__ro = pp->p.wx * pp->n * rad(0,0,1);
    pp->wy_n__ro = pp->p.wy * pp->n * rad(0,0,1);
    pp->wz__ro = pp->p.wz * rad(0,0,1);
    pp->ms_e2 = pp->p.ms * pp->e2;
}

void translateHelmertInv(struct HelmertParam *pp,
        double lat, double lon, double h, double *latp, double *lonp) {
    double sin_lat, cos_lat;
    double sin_lon, cos_lon;
    double q, n;

    if (unlikely(!pp)) {
        *latp = lat;
        *lonp = lon;
        return;
    }
    
    sin_lat = sin(lat);
    cos_lat = cos(lat);
    sin_lon = sin(lon);
    cos_lon = cos(lon);
    q = 1 / (1 - pp->e2 * sin_lat * sin_lat);
    n = pp->a * sqrt(q);

   *latp = lat
        - ((n * (q * pp->de2__2 + pp->dxe2__2) * sin_lat + pp->p.dz) * cos_lat
           - (pp->p.dx * cos_lon + pp->p.dy * sin_lon) * sin_lat
          ) / (n * q * pp->n + h)
        + (pp->wx_2e2__ro * sin_lon - pp->wy_2e2__ro * cos_lon)
          * (cos_lat * cos_lat + pp->n__2e2)
        + pp->ms_e2 * sin_lat * cos_lat;
    *lonp = lon
        + ((pp->p.dx * sin_lon - pp->p.dy * cos_lon) / (n + h)
           - (pp->wx_n__ro * cos_lon + pp->wy_n__ro * sin_lon) * sin_lat
          ) / cos_lat
        + pp->wz__ro;
}

这些都是从哪里来的？可以在proj4项目中找到一种更易于理解的语言，但是我需要优化执行速度，然后根据公式对角度总和的正弦进行任何计算，通过括号内的饰面优化指数，并分别计算常数。

现在，为了接近完成创建图块的原始任务，我们需要考虑一个称为WebMercator的坐标系。此坐标系用于OsmAnd应用程序和Web中，例如，在Google地图和OpenStreetMap中。 WebMercator是在球体上构建的Mercator投影。此投影中的坐标是像素X，Y的坐标，它们取决于Z比例，对于零比例的整个地球表面（最多约85度纬度）放置在一个256x256像素图块上，X，Y坐标从左开始从0更改为255上角，比例1-已经是4个图块，X，Y-从0到511，依此类推。

以下功能用于从WebMercator转换为WGS84：

void XYZ_WGS84(unsigned x, unsigned y, int z, double *latp, double *lonp) {
    double s = M_PI / ((1UL << 7) << z);
    *lonp = s * x - M_PI;
    *latp = asin(tanh(M_PI - s * y));
}
void WGS84_XYZ(int z, double lat, double lon, unsigned *xp, unsigned *yp) {
    double s = ((1UL << 7) << z) / M_PI;
    *xp = uint_round((lon + M_PI) * s);
    *yp = uint_round((M_PI - atanh(sin(lat))) * s);
}

在本文的第一部分结尾，我们已经可以勾勒出创建图块的算法的开始。每个256x256像素的图块使用三个值寻址：x，y，z，与坐标X，Y，Z的关系非常简单：x =（int）（X / 256）；y =（int）（Y / 256）; z = Z;

void renderTile(int z, unsigned long x, unsigned long y) {
    int i, j;
    double wlat, wlon;
    double lat, lon;

    for (i = 0; i < 255; ++i) {
        for (j = 0; j < 255; ++j) {
            XYZ_WGS84(x * 256 + j, y * 256 + i, z, &wlat, &wlon);
            translateHelmertInv(&Helmert_SK42_WGS84, wlat, wlon, 0, &lat, &lon);
            /* lat,lon -   42 */
        }
    }
}

SK42中的坐标已经转换为地图坐标系，现在仍然可以通过这些坐标在地图上找到像素，并将其颜色输入到坐标j（即i）处的图块像素中。这将是下一篇文章，其中我们将讨论测地线投影和仿射变换。

从栅格地图创建图块

More articles: