在上一部分中，我们停止了这样一个事实，即必须通过
已经转换为该地图的SC的大地坐标（纬度，经度）从原始栅格地图中获取像素的颜色。
大地坐标设置在椭球表面上，像素坐标设置在平面上，即您需要一种从凸面到平坦表面的方法。问题在于，凸面（假设已应用某种图形或坐标网格）无法无任何变形地转移到平面上。可能会变形：形状（角度），面积，线性尺寸。当然，有可能（而不是唯一一种）转移到平坦的表面上而不会仅扭曲一件事，但是其余的不可避免地会扭曲。
显然，在较小的比例（整个星球，整个大陆）上，失真比在较大的比例（区域，城市）上更明显，而在最大的比例（小区域计划）上，则完全没有讨论，因为这样的比例尺上的椭球表面不再与平面区分开。
这里我们来谈谈地图投影的概念。我将不提供将球体投影到圆柱体或圆锥体上并进行后续开发的图片示例，对我们而言，现在进行归纳非常重要。
地图投影是将椭圆形表面映射到平面上的数学定义方式。简而言之，这些是将大地坐标（纬度，经度）转换为平面笛卡尔坐标的数学公式，这正是我们所需要的。
已经发明了各种各样的制图投影，它们都为自己的目的找到了应用：用于政治地图，土壤图，等角图（在其中保留形状）的等大小（保留区域）-用于导航，等距（在选定方向上保持比例）-用于计算机处理和存储地理数据数组。也有将上述特征按一定比例组合在一起的投影，完全是深奥的投影。可以在Wikipedia的“地图投影列表”页面上找到投影示例。
对于每个投影，要么导出精确的公式，要么以无穷收敛序列的和的形式导出，在后一种情况下，甚至可能有多个选择。投影公式将纬度和经度转换为笛卡尔坐标，通常在此类坐标中以米为单位。有时可以在地图上绘制此一米长的矩形网格（以公里网格的形式），但在大多数情况下不会绘制。但是我们现在知道它仍然以隐式形式存在。只是相对于此网格的大小来计算地图上显示的地图比例尺。应当清楚地理解，这种坐标网格上的1米并不对应于地面上的1米，而并非在整个地图上，而是仅在某些点，沿着一条直线或沿着某个方向的直线上，在其余的点或方向上，会出现变形，并且地面上的1米不对应于1米的坐标网格。
一个小题外话。 WGS84_XYZ文章第一部分的功能恰恰是将WGS84 SC的坐标转换为直角坐标，但不是以米为单位，而是以像素为单位。如您所见，公式非常简单，如果墨卡托投影不是用在球体上，而是用在椭圆体上，那么公式将更加复杂。这就是为什么为了使浏览器的生活更轻松而选择了一个范围，随后WebMercator投影在其范围内扎根，对此经常受到批评。
为了我的需要，我使用了pdf格式的名为“美国地质调查局使用的地图投影”的文档，该文档可以在Internet上找到。该文档为每个投影提供了详细的说明，以简化用编程语言编写转换函数的过程。
让我们继续编写我们的算法。让我们实现一种流行的投影，称为“横轴墨卡托”，及其变体之一，称为高斯-克鲁格投影。
struct TransverseMercatorParam {
struct Ellipsoid *ep;
double k; /* */
double lat0; /* ( ) */
double lon0; /* ( ) */
double n0; /* */
double e0; /* */
double zw; /* ( ) */
double zs; /* */
//
double e2__a2k2, ie2__a2k2, m0, mp, imp;
double f0, f2, f4, f6;
double m1, m2, m4, m6;
double q, q1, q2, q3, q4, q6, q7, q8;
double q11, q12, q13, q14, q15, q16, q17, q18;
// - 2
double apk, n, b, c, d;
double b1, b2, b3, b4;
};
struct TransverseMercatorParam TM_GaussKruger = {
.ep = &Ellipsoid_Krasovsky,
.zw = rad(6,0,0),
.lon0 = -rad(3,0,0),
.e0 = 5e5,
.zs = 1e6,
};
墨卡托横向投影的一个特征是它是保形的，也就是说，地图上物体的形状和角度得以保留，并且比例尺沿一个中央子午线得以保留。该投影适用于整个地球，但是变形会随着距中央子午线的距离而增长，因此整个地球的表面会沿着子午线切成狭窄的条带，称为区域，每个子午线都使用自己的中央子午线。这种投影的实施方式的例子是高斯-克鲁格投影和UTM，其中还定义了在区域上分配坐标的方法，即：定义和相同的名称SC。
并且，实际上，是代码的初始化和坐标转换功能。在初始化函数中，常数被一次性计算，将被转换函数重用。
void setupTransverseMercator(struct TransverseMercatorParam *pp) {
double sin_lat, cos_lat, cos2_lat;
double q, n, rk, ak;
if (!pp->k)
pp->k = 1.0;
sin_lat = sin(pp->lat0);
cos_lat = cos(pp->lat0);
cos2_lat = cos_lat * cos_lat;
q = pp->ep->e2 / (1 - pp->ep->e2);
// n = (a-b)/(a+b)
n = (pp->ep->a - pp->ep->b) / (pp->ep->a + pp->ep->b);
rk = (pp->ep->a + pp->ep->b) * pp->k / 2;
ak = pp->ep->a * pp->k;
pp->e2__a2k2 = pp->ep->e2 / (ak * ak);
pp->ie2__a2k2 = (1 - pp->ep->e2) / (ak * ak);
pp->f6 = 1097.0/4 * n*n*n*n;
pp->f4 = (151.0/3 - 3291.0/8 * n) * n*n*n;
pp->f2 = (21.0/2 + (-151.0/3 + 5045.0/32 * n) * n) * n*n;
pp->f0 = (3.0 + (-21.0/4 + (31.0/4 - 657.0/64 * n) * n) * n) * n;
pp->m6 = rk * 315.0/4 * n*n*n*n;
pp->m4 = rk * (-70.0/3 - 945.0/8 * n) * n*n*n;
pp->m2 = rk * (15.0/2 + (70.0/3 + 1515.0/32 * n) * n) * n*n;
pp->m1 = rk * (-3.0 + (-15.0/4 + (-4.0 - 255.0/64 * n) * n) * n) * n;
// polar distance
pp->mp = rk * (1.0 + (1.0/4 + 1.0/64 * n*n) * n*n);
pp->imp = 1 / pp->mp;
pp->m0 = pp->n0 - pp->mp * pp->lat0 - sin_lat * cos_lat *
(pp->m1 + (pp->m2 + (pp->m4 + pp->m6 * cos2_lat) * cos2_lat) * cos2_lat);
pp->q = q;
pp->q1 = 1.0/6 * q*q;
pp->q2 = 3.0/8 * q;
pp->q3 = 5.0/6 * q;
pp->q4 = 1.0/6 - 11.0/24 * q;
pp->q6 = 1.0/6 * q;
pp->q7 = 3.0/5 * q;
pp->q8 = 1.0/5 - 29.0/60 * q;
pp->q11 = - 5.0/12 * q;
pp->q12 = -5.0/24 + 3.0/8 * q;
pp->q13 = - 1.0/240 * q*q;
pp->q14 = 149.0/360 * q;
pp->q15 = 61.0/720 - 63.0/180 * q;
pp->q16 = - 1.0/40 * q*q;
pp->q17 = - 1.0/60 * q;
pp->q18 = 1.0/24 + 1.0/15 * q;
// - 2
double e2 = pp->ep->e2;
pp->apk = ak * (1 + n*n / 4 + n*n*n*n / 64) / (1 + n);
pp->n = n;
pp->b = (5 - e2) * e2 * e2 / 6;
pp->c = (104 - 45 * e2) * e2 * e2 * e2 / 120;
pp->d = 1237.0/1260 * e2 * e2 * e2 * e2;
pp->b1 = (1.0/2 + (-2.0/3 + (5.0/16 + 41.0/180 * n) * n) * n) * n;
pp->b2 = (13.0/48 + (-3.0/5 + 557.0/1440 * n) * n) * n*n;
pp->b3 = (61.0/240 - 103.0/140 * n) * n*n*n;
pp->b3 = 49561.0/161280 * n*n*n*n;
}
void translateTransverseMercator(struct TransverseMercatorParam *pp, int zone,
double lat, double lon, double *ep, double *np) {
double lon2, v, m;
double k4, k6, h3, h5;
double sin_lat = sin(lat);
double cos_lat = cos(lat);
double cos2_lat = cos_lat * cos_lat;
lon -= zone * pp->zw + pp->lon0;
while (unlikely(lon <= -M_PI))
lon += 2*M_PI;
lon2 = lon * lon;
//
v = 1 / sqrt(pp->e2__a2k2 * cos2_lat + pp->ie2__a2k2);
m = ((pp->m6 * cos2_lat + pp->m4) * cos2_lat + pp->m2) * cos2_lat + pp->m1;
k4 = ((pp->q1 * cos2_lat + pp->q2) * cos2_lat + 1.0/4 ) * cos2_lat - 1.0/24;
k6 = ((pp->q3 * cos2_lat + pp->q4) * cos2_lat - 1.0/12) * cos2_lat + 1.0/720;
h3 = (( pp->q6) * cos2_lat + 1.0/3 ) * cos2_lat - 1.0/6;
h5 = ((pp->q7 * cos2_lat + pp->q8) * cos2_lat - 1.0/6 ) * cos2_lat + 1.0/120;
// ( )
*np = pp->m0 + pp->mp * lat
+ (m + v * ((k6 * lon2 + k4) * lon2 + 0.5) * lon2) * cos_lat * sin_lat;
*ep = pp->e0 + pp->zs * zone
+ ( v * ((h5 * lon2 + h3) * lon2 + 1.0) * lon ) * cos_lat;
}
在变换函数的输出处，我们将有坐标：东和北位移（e，n）是以米为单位的直角笛卡尔坐标。
我们已经非常接近完成我们的算法。我们只需要在地图文件中找到像素（x，y）的坐标即可。因为 像素坐标也是笛卡尔坐标，那么我们可以通过仿射变换（e，n）到（x，y）来找到它们。稍后，我们将返回查找最仿射变换的参数。
struct AffineParam {
double c00, c01, c02;
double c10, c11, c12;
};
void translateAffine(struct AffineParam *app, double e, double n,
double *xp, double *yp) {
*xp = app->c00 + app->c01 * e + app->c02 * n;
*yp = app->c10 + app->c11 * e + app->c12 * n;
}
最后，完整的图块创建算法：
void renderTile(ImagePtr tile, int z, unsigned long x, unsigned long y) {
int i, j;
double wlat, wlon;
ImagePtr srcimg;
double lat, lon;
double e, n;
double x, y;
for (i = 0; i < 256; ++i) {
for (j = 0; j < 256; ++j) {
XYZ_WGS84(x * 256 + j, y * 256 + i, z, &wlat, &wlon);
translateHelmertInv(&Helmert_SK42_WGS84, wlat, wlon, 0, &lat, &lon);
findSrcImg(&srcimg, lat, lon);
translateTransverseMercator(&TM_GaussKruger, srcimg->zone, lat, lon, &e, &n);
translateAffine(&srcimg->affine, e, n, &x, &y);
setPixelColor(tile, j, i, interpolatePixelColor(srcimg, x, y));
}
}
}
z = 12，y = 1377，x = 2391的工作结果：
在该算法中，仍未描述从地图的SC中指定的大地坐标lat，lon查找srcimg地图的原始图像的功能。我认为它和srcimg-> zone zone的数量没有问题，但是我们将详细介绍仿射变换srcimg->仿射的参数。
很久以前，在Internet上的某个地方，发现了一个函数来查找仿射变换的参数，即使引用原始注释，我也引用它：
struct TiePoint {
struct TiePoint *next;
double lon, lat;
double e, n;
double x, y;
};
void setupAffine(struct AffineParam *app, struct TiePoint *plist) {
/*
* :
* x = c00 + c01 * e + c02 * n
* y = c10 + c11 * e + c12 * n
*/
struct TiePoint *pp; /* */
double a11, a12, a13; /* */
double a21, a22, a23; /* 3*3 */
double a31, a32, a33; /* */
double b1, b2, b3; /* */
int m; /* z: m=0 -> z=x, m=1 -> z=y */
double z; /* x, y */
double t; /* */
/* 2- 3 ,
. */
/* */
a11 = a12 = a13 = a22 = a23 = a33 = 0;
for (pp = plist; pp; pp = pp->next) {
a11 += 1;
a12 += pp->e;
a13 += pp->n;
a22 += pp->e * pp->e;
a23 += pp->e * pp->n;
a33 += pp->n * pp->n;
}
/* ( ) */
a21 = a12;
a31 = a13;
a12 /= a11;
a13 /= a11;
a22 -= a21 * a12;
a32 = a23 - a31 * a12;
a23 = a32 / a22;
a33 -= a31 * a13 + a32 * a23;
/* , X Y */
for (m = 0; m < 2; m++) { /* m=0 -> X, m=1 -> Y */
/* */
b1 = b2 = b3 = 0;
for (pp = plist; pp; pp = pp->next) {
z = !m ? pp->x : pp->y;
b1 += z;
b2 += pp->e * z;
b3 += pp->n * z;
}
/* */
b1 /= a11;
b2 = (b2 - a21 * b1) / a22;
b3 = (b3 - a31 * b1 - a32 * b2) / a33;
/* */
t = b2 - a23 * b3;
*(!m ? &app->c00 : &app->c10) = b1 - a12 * t - a13 * b3;
*(!m ? &app->c01 : &app->c11) = t;
*(!m ? &app->c02 : &app->c12) = b3;
}
}
在输入处，至少必须提交三个锚点，在输出处，我们获得现成的参数。锚点是同时已知像素坐标（x，y）和东，北偏移坐标（e，n）的点。原始地图上公里网格的交点可以用作此类点。如果地图上没有公里网格怎么办？然后，您可以设置对（x，y）和（lon，lat），因为这些点采用了平行线和子午线的交点，因此它们始终在地图上。它仅保留将（lon，lat）转换为（e，n），这是通过投影的转换函数完成的，在我们的例子中是translateTransverseMercator（）。
如您所见，需要锚点来告诉算法带有地图图像的文件描述了地球表面的哪一部分。由于两个坐标系都是笛卡尔坐标系，因此无论我们设置多少个锚点，以及它们彼此之间的距离如何，整个地图平面上的差异都将在确定锚点的误差之内。大部分错误是由于使用了错误的（没有指定参数）投影，基准面或椭球，因此在笛卡尔坐标系中无法获得输出的坐标（e，n），而是相对于笛卡尔坐标系略有弯曲。在视觉上，这可以看成是“碎纸”。自然地，增加锚点的数量并不能解决这个问题。可以通过调整投影，基准面和椭球面的参数来解决该问题，在这种情况下，大量的锚点将使您更准确地平滑“工作表”，而不会错过未平滑的区域。
在文章的结尾，我将告诉您如何在OpenLayers中使用现成的图块，以及以何种形式为OsmAnd程序做准备。
对于OpenLayers，您只需要将磁贴放置在Web上并命名即可，以便可以在文件或目录名称中突出显示（z，y，x），例如：
/tiles/topo/12_1377_2391.jpg
或什至更好：
/ tile / topo /12/1377/2391.jpg
然后可以像这样使用它们：
map = new OpenLayers.Map (...);
map.addLayer(new OpenLayers.Layer.XYZ("Topo Maps", "/tiles/topo/${z}_${y}_${x}.jpg", {
isBaseLayer: false, visibility: false,
}));
对于OsmAnd程序，很容易从具有一组磁贴的任何现有文件中确定格式。我马上告诉你结果。这些磁贴被打包到如下创建的sqlite数据库文件中：
CREATE TABLE info AS SELECT 99 AS minzoom, 0 AS maxzoom;
CREATE TABLE tiles (x int, y int, z int, s int, image blob, PRIMARY KEY (x,y,z,s));
CREATE INDEX IND on tiles (x,y,z,s);
UPDATE info SET minzoom=$minzoom, maxzoom=$maxzoom
列s总是用零填充，对此，我不明白，图像以原始二进制形式输入，格式（jpg，png，gif）丢失了，但是OsmAnd由其内容确定。可以将不同格式的图块打包到一个数据库中。代替$ minzoom和$ maxzoom，有必要替换输入到基础中的图块的比例限制。已完成的数据库文件（例如Topo.sqlitedb）被转移到osmand / tiles目录中的智能手机或平板电脑。重新启动OsmAnd，在菜单->“配置地图”->“顶层”中，将出现一个新选项“ Topo”-这是带有新图块的地图。