前の部分では
、このマップのSCにすでに変換されている地理座標（緯度、経度）によって、元のラスターマップからピクセルの色を取得する必要があるという事実にとどまりました。
ジオデティック座標は楕円形の表面に設定され、ピクセル座標は平面上にあります。凸面から平らな面に移行する方法が必要です。問題は、凸面（ある種の描画または座標グリッドが適用されていると想像してください）を歪みなしに平面に転写できないことです。歪む可能性があります：形状（角度）、面積、直線寸法。もちろん、1つだけを歪めることなく平らな面に移動する可能性はありますが、残りは必然的に歪められます。
明らかに、小規模（惑星全体、大陸）では、歪みは大規模なもの（地域、都市）よりも顕著であり、最大規模（小さな領域の計画）では、歪みについてはまったく話していません。そのようなスケールの楕円形の表面は、もはや平面と区別できません。
ここで、マップ投影の概念について説明します。球を円柱または円錐に投影し、その後に展開する写真の例は示しません。ここで一般化することが重要です。
マッププロジェクションは、楕円形の表面を平面にマッピングする数学的に定義された方法です。簡単に言えば、これらは測地座標（緯度、経度）を平坦なカルテシアン座標に変換するための数式です-まさに必要なものです。
多種多様なカートグラフィックプロジェクションが発明され、それらはすべて独自の目的に適用されます：政治地図、土壌地図、コンフォーマル（形状が保存される）-ナビゲーション用、等距離（選択された方向にスケールを維持する）-コンピュータ用の同じサイズ（領域が保存される）ジオデータ配列の処理と保存。上記の特徴を一定の比率で組み合わせたプロジェクションもあり、完全に難解なプロジェクションがあります。投影法の例は、Wikipediaの「マップ投影法のリスト」ページにあります。
射影ごとに、正確な式が導き出されるか、無限の収束系列の合計の形で導き出されます。後者の場合、いくつかのオプションがあります。射影式は緯度と経度をカルテシアン座標に変換します。通常、メーターはそのような座標の測定単位として使用されます。この1メートルの長方形グリッドは、マップ上に（キロメートルグリッドの形式で）プロットできる場合がありますが、ほとんどの場合、プロットされません。しかし、それがまだ暗黙の形で存在していることがわかりました。マップに示されているマップの縮尺は、このグリッドのサイズを基準にして計算されます。このような座標グリッド上の1メートルは、マップ全体ではなく、特定の線に沿った、または特定の方向の線に沿った特定のポイントでのみ、地上の1メートルに対応することを明確に理解する必要があります。他のポイントまたは方向では歪みが発生し、地面の1メートルは座標グリッドの1メートルに対応しません。
小さな逸脱。記事WGS84_XYZの最初の部分の関数は、正確にはWGS84 SCから長方形の座標への座標の変換ですが、メートルではなくピクセルで表されます。ご覧のとおり、メルカトル投影が球ではなく楕円形で使用された場合、式は非常に単純になり、式はより複雑になります。そのため、ブラウザの作業を楽にするために球が選択され、その後WebMercatorプロジェクションがその球に根付いており、しばしば批判されています。
必要に応じて、インターネットで入手できる「米国地質調査で使用された地図投影法」というpdf形式のドキュメントを使用しました。このドキュメントでは、プログラミング言語で変換関数を簡単に記述できるように、各プロジェクションの詳細な手順を説明しています。
アルゴリズムを書き続けましょう。 Transverse Mercatorと呼ばれる人気のあるプロジェクションの1つと、Gauss-Krugerプロジェクションと呼ばれるそのバリアントの1つを実装しましょう。
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,
};
横方向のメルカトル投影の特徴は、それがコンフォーマルであるということです。つまり、マップ上のオブジェクトの形状と角度が保持され、スケールが1つの中央子午線に沿って保持されるという事実です。投影は地球全体に適していますが、中央子午線からの距離とともに歪みが大きくなるため、地表全体が子午線に沿って狭い帯状に切断されます。これはゾーンと呼ばれ、それぞれに独自の中央子午線が使用されます。このようなプロジェクションの実装例は、Gauss-KrugerプロジェクションとUTMであり、ゾーン全体の座標の分布方法も定義されています。定義され、同じ名前のSC。
そして、実際には、座標の初期化と変換の関数のコード。初期化関数では、定数が1回計算され、変換関数によって再利用されます。
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->ゾーンゾーンの数に問題はないと思いますが、アフィン変換srcimg-> affineのパラメーターをより詳細に見つけることに専念します。
非常に昔のインターネットのどこかで、そのような関数がアフィン変換のパラメーターを見つけることがわかったので、元のコメントでも引用します。
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;
}
}
入力では、少なくとも3つのアンカーポイントを送信する必要があります。出力では、既製のパラメーターを取得します。アンカーポイントは、ピクセル座標（x、y）と東および北のオフセット座標（e、n）の両方がわかっているポイントです。元のマップ上のキロメートルグリッドの交点をそのような点として使用できます。マップ上にキロメートルグリッドがない場合はどうなりますか？次に、ペア（x、y）と（lon、lat）を設定できます。このようなポイントは、平行線と子午線の交点を取り、常にマップ上にあります。 （lon、lat）を（e、n）に変換するだけです。これは、射影の変換関数によって行われます。この場合はtranslateTransverseMercator（）です。
ご覧のとおり、アンカーポイントは、マップイメージを含むファイルが地球の表面のどの部分を記述しているかをアルゴリズムに伝えるために必要です。両方の座標系はデカルトであるため、設定したアンカーポイントの数や距離に関係なく、マッププレーン全体の不一致は、アンカーポイントを決定する際のエラーの範囲内になります。エラーのほとんどは、間違った（正確に指定されていないパラメーターを使用した）投影、データム、または楕円形が使用されていることです。その結果、出力の座標（e、n）は、カルテシアン座標系ではなく、カルテシアン座標系に対してわずかに湾曲して取得されます。視覚的には、これは「しわくちゃのシート」として視覚化できます。当然、アンカーポイントの数を増やしてもこの問題は解決しません。この問題は、投影、データム、楕円形のパラメータを調整することで解決できます。この場合、多数のアンカーポイントを使用すると、「シート」をより正確に滑らかにし、滑らかでない領域を見逃さないようにすることができます。
そして、記事の最後に、OpenLayersで既製のタイルを使用する方法と、OsmAndプログラム用にそれらを準備するための形式について説明します。
OpenLayersの場合、ファイルまたはディレクトリ名で（z、y、x）を強調表示できるように、タイルをWebに配置して名前を付ける必要があります。例：
/tiles/topo/12_1377_2391.jpg
またはさらに良い方法：
/ tiles / 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はその内容によって決定します。さまざまな形式のタイルを1つのデータベースにパックできます。$ minzoomと$ maxzoomの代わりに、ベースに入力されたタイルのスケール制限を置き換える必要があります。完成したデータベースファイル（Topo.sqlitedbなど）は、osmand / tilesディレクトリのスマートフォンまたはタブレットに転送されます。OsmAndを再起動し、[メニュー]-> [マップの構成]-> [トップレイヤー]に新しいオプション[Topo]が表示されます。これは新しいタイルを含むマップです。