この記事では、私の「研究」の結果について説明し、任意のラスターマップをアプリケーションが理解できるタイルに変換するためのアルゴリズムを作成し、その過程で、楕円形、データム、座標系、投影などの概念に精通します。
私たちの地球は球の形をしておらず、楕円形でさえありません。この複雑な図はジオイドと呼ばれています。事実、地球内部の質量は均等に分布していないため、地球がわずかに凹んでいる場所もあれば、わずかに凸状になっている場所もあります。別の国または大陸の領土を取得する場合、この楕円形の中心が地球の質量中心に対して3つの座標軸に沿ってわずかにシフトしていれば、十分な精度でその表面を楕円形の表面と位置合わせできます。このような楕円形は参照楕円形と呼ばれ、それが作成された地球のローカル領域のみを記述するのに適しています。このサイトから遠く離れていると、地球の表面と非常に大きな差異が生じる可能性があります。中心が地球の質量中心と一致する楕円形は、一般的な地上楕円形と呼ばれます。晴れ、参照楕円形は、一般的な地上よりも地球のその局所領域をよりよく説明していますが、一般的な地上は地球の表面全体に適しています。
楕円形を説明するには、赤道半径（通常はaで示されます）と極半径（b）の2つの独立した値だけで十分ですが、2番目の独立した値の代わりに、極収縮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,
};
この例は、2つの楕円形を示しています。衛星ナビゲーションで使用される一般的な地上WGS84と、ヨーロッパとアジアに適用可能なKrasovsky参照楕円形です。
もう一つの興味深い点を考えてみてください。地球の形は遅いですが変化しているので、今日表面をよく表している楕円形は、100年後には現実からかけ離れているかもしれません。これは、一般的な地上の楕円形とはほとんど関係がありません。同じエラー内の偏差ですが、参照楕円に関連しています。ここで、別の概念、つまりデータムについて説明します。データムは、特定の時点で固定された、楕円（a、f）、地球内部の変位（参照楕円の場合）のパラメーターのセットです。より正確には、データムは必ずしも楕円形を表すとは限らず、より複雑な図形、たとえば準ジオイドのパラメータである可能性があります。
確かに、ある楕円形またはデータムから別の楕円形またはデータムに移動する方法については、すでに疑問が生じています。これを行うには、各楕円形に測地座標系（緯度と経度（ファイ、ラムダ））が必要です。遷移は、座標をある座標系から別の座標系に変換することによって実行されます。
座標を変換するためのさまざまな式があります。最初に、ある座標系の測地座標を3次元座標X、Y、Zに変換し、それらを使用してシフトとターンを実行してから、結果の3次元座標を別の座標系の測地座標に変換できます。あなたはそれを直接行うことができます。なぜならすべての式は無限に収束するシリーズであり、必要な精度を達成するために、通常はシリーズのいくつかのメンバーに制限されます。例として、Helmert変換を使用します。これらは、3次元座標への遷移を使用した変換であり、上記の3つの段階で構成されます。変換には、2つの楕円形に加えて、さらに7つのパラメーターが必要です。3つの軸に沿った3つのシフト、3つの回転角度、およびスケール係数です。ご想像のとおり、すべてのパラメーターはデータムから抽出できます。ただし、アルゴリズムでは、このようなオブジェクトをデータムとして使用せず、代わりに、ある座標系から別の座標系への遷移オブジェクトを導入します。これには、2つの楕円体への参照と7つの変換パラメーターが含まれます。これは、アルゴリズムの2番目のオブジェクトになります。
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},
};
この例は、Pulkovo1942座標系からWGS84座標系に変換するためのパラメーターを示しています。変換パラメーター自体は別のトピックです。開いている座標系もあれば、経験的に選択された座標系もあるため、ソースによって値がわずかに異なる場合があります。
パラメータに加えて、変換関数も必要です。これは直接であり、反対方向への変換の場合は、反対方向への変換のみが必要です。たくさんの数学をスキップして、最適化された関数を提供します。
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度の緯度）が1つの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ピクセルの各タイルは、次の3つの値でアドレス指定されます。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のタイルピクセルに入力する必要があります。これは次の記事で、ジオデシックプロジェクションとアフィン変換について説明します。