Dalam artikel ini saya akan memberi tahu Anda tentang hasil "penelitian" saya, menyusun algoritme untuk mengubah peta raster arbitrer menjadi ubin yang dapat dipahami untuk aplikasi dan, di sepanjang jalan, berkenalan dengan konsep seperti elipsoid, datum, sistem koordinat, proyeksi.
Bumi kita tidak berbentuk bola, dan bahkan tidak berbentuk elipsoid; Sosok kompleks ini disebut geoid. Faktanya adalah bahwa massa di dalam bumi tidak terdistribusi secara merata, sehingga di beberapa tempat bumi agak cekung, di tempat lain agak cembung. Jika kita mengambil wilayah negara atau benua yang terpisah, maka permukaannya dengan ketelitian yang memadai dapat disejajarkan dengan permukaan beberapa ellipsoid, jika pusat ellipsoid ini sedikit bergeser sepanjang tiga sumbu koordinat relatif terhadap pusat massa bumi. Elipsoid semacam itu disebut ellipsoid referensi, sangat cocok untuk mendeskripsikan hanya area lokal Bumi tempat ia dibuat. Pada jarak yang jauh dari situs ini, ia dapat memiliki perbedaan yang sangat besar dengan permukaan bumi. Elipsoid yang pusatnya bertepatan dengan pusat massa Bumi disebut ellipsoid terestrial umum. Bersih,bahwa ellipsoid referensi lebih baik menggambarkan bagian lokalnya di Bumi daripada terestrial umum, tetapi terestrial umum cocok untuk seluruh permukaan Bumi.
Untuk mendeskripsikan elipsoid, hanya dua nilai independen yang cukup: jari-jari ekuator (biasanya dilambangkan dengan a) dan jari-jari kutub (b), tetapi alih-alih nilai independen kedua, kontraksi kutub f = (ab) / a biasanya digunakan. Ini adalah hal pertama yang kita butuhkan dalam algoritme kita sebagai objek - elipsoid. Untuk bagian bumi yang berbeda pada tahun yang berbeda, peneliti yang berbeda telah menghitung banyak elipsoid referensi, informasi tentangnya diberikan dalam bentuk data: a (dalam meter) dan 1 / f (tanpa dimensi). Anehnya, untuk ellipsoid terestrial biasa, ada juga banyak varian yang berbeda (berbeda a, f), tetapi perbedaannya tidak terlalu kuat, ini terutama disebabkan oleh perbedaan metode untuk menentukan a dan 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,
};
Contoh ini menunjukkan dua ellipsoid: WGS84 terestrial umum, digunakan dalam navigasi satelit, dan ellipsoid referensi Krasovsky, berlaku untuk Eropa dan Asia.
Pertimbangkan hal menarik lainnya, faktanya adalah bahwa bentuk bumi lambat, tetapi berubah, sehingga ellipsoid, yang sekarang menggambarkan permukaan dengan baik, dalam seratus tahun mungkin jauh dari kenyataan. Ini tidak ada hubungannya dengan ellipsoid terestrial biasa, sejak itu deviasi dalam kesalahan yang sama, tetapi relevan untuk referensi ellipsoid. Di sini kita sampai pada konsep lain - datum. Datum adalah sekumpulan parameter ellipsoid (a, f), perpindahannya di dalam Bumi (untuk referensi ellipsoid), ditetapkan pada titik waktu tertentu. Lebih tepatnya, datum belum tentu menggambarkan elipsoid, itu bisa menjadi parameter dari sosok yang lebih kompleks, misalnya, kuasigeoid.
Tentunya pertanyaan sudah muncul: bagaimana cara berpindah dari satu elipsoid atau datum ke yang lain? Untuk ini, setiap elipsoid harus memiliki sistem koordinat geodetik: lintang dan bujur (phi, lambda), transisi dilakukan dengan menerjemahkan koordinat dari satu sistem koordinat ke sistem lainnya.
Ada berbagai rumus untuk mengubah koordinat. Pertama-tama Anda dapat menerjemahkan koordinat geodesik dalam satu sistem koordinat menjadi koordinat tiga dimensi X, Y, Z, melakukan pergeseran dan rotasi dengannya, lalu mengubah koordinat tiga dimensi yang dihasilkan menjadi koordinat geodetik di sistem koordinat lain. Anda bisa melakukannya secara langsung. Karena Karena semua rumus adalah deret konvergen tak hingga, biasanya dibatasi pada beberapa anggota deret untuk mencapai akurasi yang diperlukan. Sebagai contoh, kita akan menggunakan transformasi Helmert, yaitu transformasi menggunakan transisi ke koordinat tiga dimensi, terdiri dari tiga tahap yang dijelaskan di atas. Untuk transformasi, selain dua elipsoid, diperlukan 7 parameter lagi: tiga pergeseran sepanjang tiga sumbu, tiga sudut rotasi, dan faktor skala. Seperti yang Anda duga, semua parameter dapat diekstraksi dari datum.Namun dalam algoritme kami tidak akan menggunakan objek seperti datum, tetapi kami akan memperkenalkan objek transisi dari satu sistem koordinat ke sistem koordinat lainnya, yang akan berisi: referensi ke dua elipsoid dan 7 parameter transformasi. Ini akan menjadi objek kedua dari algoritme kami.
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},
};
Contoh ini menunjukkan parameter untuk mengubah sistem koordinat Pulkovo 1942 menjadi sistem koordinat WGS84. Parameter transformasi itu sendiri adalah topik yang terpisah, untuk beberapa sistem koordinat terbuka, untuk yang lain mereka dipilih secara empiris, sehingga nilainya mungkin sedikit berbeda di sumber yang berbeda.
Selain parameter, diperlukan juga fungsi transformasi, bisa langsung dan untuk transformasi ke arah berlawanan, kita hanya perlu transformasi ke arah yang berlawanan. Saya akan melewatkan banyak matematika dan memberikan fungsi saya yang dioptimalkan.
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;
}
Darimana semua ini berasal? Dalam bahasa yang lebih dimengerti, rumus dapat ditemukan di proyek proj4, tetapi sejak itu Saya membutuhkan pengoptimalan untuk kecepatan eksekusi, kemudian setiap perhitungan sinus dari jumlah sudut diubah oleh rumus, eksponen dioptimalkan oleh veneing dalam tanda kurung, dan konstanta dihitung secara terpisah.
Sekarang, untuk lebih dekat dalam menyelesaikan tugas awal membuat ubin, kita perlu mempertimbangkan sistem koordinat yang disebut WebMercator. Sistem koordinat ini digunakan di aplikasi OsmAnd dan di web, misalnya, di Google maps dan di OpenStreetMap. WebMercator adalah proyeksi Mercator yang dibangun di atas bola. Koordinat dalam proyeksi ini adalah koordinat piksel X, Y dan bergantung pada skala Z, untuk skala nol seluruh permukaan bumi (hingga sekitar 85 derajat lintang) ditempatkan pada satu petak piksel 256x256, koordinat X, Y berubah dari 0 menjadi 255, dimulai dari kiri pojok atas, untuk skala 1 - sudah 4 ubin, X, Y - dari 0 hingga 511 dan seterusnya.
Fungsi berikut digunakan untuk mengonversi dari WebMercator ke 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);
}
Dan di akhir bagian pertama artikel, kita sudah bisa membuat sketsa awal algoritma kita untuk membuat ubin. Setiap ubin 256x256 piksel ditangani oleh tiga nilai: x, y, z, hubungan dengan koordinat X, Y, Z sangat sederhana: 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 */
}
}
}
Koordinat di SK42 sudah diubah koordinatnya menjadi sistem koordinat peta, sekarang tinggal mencari piksel pada peta dengan koordinat tersebut dan memasukkan warnanya menjadi piksel petak pada koordinat j, i. Ini akan menjadi artikel berikutnya, di mana kita akan berbicara tentang proyeksi geodesik dan transformasi affine.