Pada bagian sebelumnya, kami berhenti pada fakta bahwa perlu mendapatkan warna piksel dari peta raster asli dengan koordinat geodetik (lintang,
bujur), yang sudah diterjemahkan ke dalam SC peta ini.
Koordinat geodetik ditetapkan pada permukaan elipsoid, dan koordinat piksel pada bidang, mis. Anda membutuhkan cara untuk beralih dari permukaan cembung ke permukaan datar. Masalahnya adalah permukaan cembung (bayangkan bahwa beberapa jenis gambar atau kisi koordinat diterapkan padanya) tidak dapat dipindahkan ke permukaan datar tanpa distorsi. Mungkin terdistorsi: bentuk (sudut), luas, dimensi linier. Tentu saja ada kemungkinan, dan bukan satu-satunya, untuk berpindah ke permukaan datar tanpa hanya mendistorsi satu hal, tetapi yang lainnya pasti akan terdistorsi.
Jelas, pada skala yang lebih kecil (seluruh planet, benua), distorsi lebih terlihat daripada pada yang lebih besar (wilayah, kota), dan pada yang terbesar (rencana area kecil), kita tidak membicarakannya sama sekali, karena permukaan ellipsoid pada skala seperti itu tidak lagi dapat dibedakan dari bidangnya.
Di sini kita sampai pada konsep proyeksi peta. Saya tidak akan memberikan contoh dengan gambar proyeksi bola ke silinder atau kerucut dengan pengembangan selanjutnya, bagi kami sekarang penting untuk menggeneralisasi.
Proyeksi peta adalah cara yang ditentukan secara matematis untuk memetakan permukaan elipsoid ke bidang. Sederhananya, ini adalah rumus matematika untuk mengubah koordinat geodetik (lintang, bujur) menjadi koordinat Kartesius datar - sesuai yang kita butuhkan.
Berbagai macam proyeksi kartografi telah ditemukan, semuanya menemukan aplikasi untuk tujuan mereka sendiri: berukuran sama (di mana area tersebut dipertahankan) untuk peta politik, peta tanah, konformal (di mana bentuknya diawetkan) - untuk navigasi, jarak yang sama (mempertahankan skala dalam arah yang dipilih) - untuk komputer pemrosesan dan penyimpanan array geodata. Ada juga proyeksi yang menggabungkan fitur-fitur yang disebutkan di atas dalam proporsi tertentu, ada juga proyeksi yang sepenuhnya esoterik. Contoh proyeksi dapat ditemukan di Wikipedia pada halaman Daftar Proyeksi Peta.
Untuk setiap proyeksi, baik rumus yang tepat diturunkan, atau dalam bentuk jumlah deret konvergen tak hingga, dan dalam kasus terakhir bahkan mungkin ada beberapa opsi. Rumus proyeksi mengubah lintang dan bujur menjadi koordinat Kartesius, biasanya meteran digunakan sebagai satuan pengukuran dalam koordinat tersebut. Bingkai persegi 1 meter seperti itu terkadang dapat diplot di peta (dalam bentuk kisi kilometer), tetapi dalam banyak kasus, bingkai tidak diplot. Tapi sekarang kita tahu bahwa itu masih ada dalam bentuk implisit. Skala peta, yang ditunjukkan pada peta, hanya dihitung relatif terhadap ukuran grid ini. Harus dipahami dengan jelas bahwa 1 meter pada kisi koordinat seperti itu sama dengan 1 meter di permukaan tanah, bukan pada keseluruhan peta, tetapi hanya pada beberapa titik, sepanjang garis tertentu, atau sepanjang garis dalam arah tertentu,di titik atau arah lain muncul distorsi dan 1 meter di permukaan tanah tidak sesuai dengan 1 meter garis koordinat.
Penyimpangan kecil. Fungsi dari bagian pertama artikel WGS84_XYZ tepatnya adalah transformasi koordinat dari WGS84 SC menjadi koordinat persegi panjang, tetapi dinyatakan bukan dalam meter, tetapi dalam piksel. Seperti yang Anda lihat, rumus di sana sangat sederhana, jika proyeksi Mercator tidak digunakan pada bola, tetapi pada elipsoid, rumusnya akan menjadi lebih rumit. Itulah sebabnya, untuk membuat hidup lebih mudah bagi browser, sebuah bola dipilih, kemudian proyeksi WebMercator berakar dengan bulatannya, yang sering dikritik.
Untuk kebutuhan saya, saya menggunakan dokumen yang disebut "Proyeksi peta yang digunakan oleh US Geological Survey" dalam format pdf, yang dapat ditemukan di Internet. Dokumen tersebut memberikan instruksi rinci untuk setiap proyeksi untuk memudahkan penulisan fungsi transformasi dalam bahasa pemrograman.
Mari lanjutkan menulis algoritma kita. Mari terapkan salah satu proyeksi populer yang disebut Transverse Mercator dan salah satu variannya yang disebut proyeksi Gauss-Kruger.
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,
};
Fitur proyeksi Mercator transversal adalah bahwa ia konformal, yaitu, bentuk objek pada peta dan sudut dipertahankan, serta fakta bahwa skala dipertahankan di sepanjang satu meridian pusat. Proyeksi ini cocok untuk seluruh globe, tetapi distorsi tumbuh dengan jarak dari meridian pusat, sehingga seluruh permukaan bumi dipotong menjadi strip sempit di sepanjang meridian, yang disebut zona, yang masing-masing menggunakan meridian tengahnya sendiri. Contoh implementasi proyeksi tersebut adalah proyeksi Gauss-Kruger dan UTM, di mana metode distribusi koordinat di atas zona juga ditentukan, yaitu. didefinisikan dan nama yang sama SC.
Dan, sebenarnya, kode untuk inisialisasi dan fungsi transformasi koordinat. Dalam fungsi inisialisasi, konstanta dihitung satu kali, yang akan digunakan kembali oleh fungsi konversi.
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;
}
Pada keluaran dari fungsi transformasi, kita akan mendapatkan koordinat: perpindahan timur dan utara (e, n) adalah koordinat Kartesius persegi panjang dalam meter.
Kami sudah sangat dekat untuk menyelesaikan algoritme kami. Kami hanya perlu menemukan koordinat piksel (x, y) di file peta. Karena koordinat piksel juga Cartesian, maka kita dapat menemukannya dengan transformasi affine (e, n) menjadi (x, y). Kami akan kembali untuk menemukan parameter transformasi affine paling banyak nanti.
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;
}
Dan terakhir, algoritma pembuatan ubin lengkap:
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));
}
}
}
Hasil pekerjaan untuk z = 12, y = 1377, x = 2391:
Dalam algoritma, fungsi untuk menemukan citra asli peta srcimg dari koordinat geodetik lat, lon yang ditentukan dalam SC peta tidak dijelaskan. Saya pikir tidak akan ada masalah dengan itu dan jumlah zona zona srcimg->, tapi kita akan memikirkan menemukan parameter dari transformasi affine srcimg-> affine secara lebih rinci.
Di suatu tempat di Internet, sangat lama sekali, fungsi seperti itu ditemukan untuk menemukan parameter transformasi affine, saya mengutipnya bahkan dengan komentar asli:
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;
}
}
Pada bagian input minimal tiga anchor point harus diserahkan, pada output kita mendapatkan parameter yang sudah jadi. Anchor point adalah titik dimana koordinat piksel (x, y) dan koordinat offset timur dan utara (e, n) diketahui. Titik-titik perpotongan dari grid kilometer pada peta asli dapat digunakan sebagai titik-titik tersebut. Bagaimana jika tidak ada kisi kilometer di peta? Kemudian Anda dapat mengatur pasangan (x, y) dan (lon, lat), karena titik-titik tersebut mengambil titik potong dari paralel dan meridian, mereka selalu berada di peta. Tetap hanya untuk mengkonversi (lon, lat) menjadi (e, n), ini dilakukan oleh fungsi transformasi untuk proyeksi, dalam kasus kami itu adalah translateTransverseMercator ().
Seperti yang Anda lihat, titik jangkar diperlukan untuk memberi tahu algoritme bagian mana dari permukaan bumi yang dijelaskan oleh file dengan gambar peta. Karena kedua sistem koordinat adalah Cartesian, tidak peduli berapa banyak titik jangkar yang kita tetapkan dan tidak peduli seberapa jauh mereka dari satu sama lain, perbedaan pada seluruh bidang peta akan menjadi kesalahan dalam menentukan titik jangkar. Sebagian besar kesalahan adalah bahwa proyeksi, datum atau ellipsoid yang salah (dengan parameter yang tidak ditentukan secara tepat) digunakan, akibatnya, koordinat (e, n) pada keluaran tidak berada dalam sistem koordinat Kartesius, tetapi dalam relatif sedikit melengkung dibandingkan dengan sistem Kartesius. Secara visual, ini bisa divisualisasikan sebagai “lembaran kusut”. Secara alami, menambah jumlah titik jangkar tidak menyelesaikan masalah ini. Masalah tersebut dapat diselesaikan dengan menyetel parameter proyeksi, datum dan ellipsoid,dalam hal ini, sejumlah besar titik jangkar akan memungkinkan Anda untuk menghaluskan "lembaran" dengan lebih akurat dan tidak melewatkan area yang tidak mulus.
Dan di akhir artikel saya akan memberi tahu Anda cara menggunakan ubin yang sudah jadi di OpenLayers dan dalam bentuk apa mempersiapkannya untuk program OsmAnd.
Untuk OpenLayers, Anda hanya perlu meletakkan ubin di web dan
menamainya sehingga Anda dapat menyorot (z, y, x) di file atau nama direktori, misalnya: /tiles/topo/12_1377_2391.jpg
atau, lebih baik lagi:
/ tiles / topo /12/1377/2391.jpg
Kemudian mereka dapat digunakan seperti ini:
map = new OpenLayers.Map (...);
map.addLayer(new OpenLayers.Layer.XYZ("Topo Maps", "/tiles/topo/${z}_${y}_${x}.jpg", {
isBaseLayer: false, visibility: false,
}));
Untuk program OsmAnd, mudah untuk menentukan format dari file yang sudah ada dengan satu set ubin. Saya akan segera memberi tahu Anda tentang hasilnya. Ubin tersebut dikemas ke dalam file database sqlite yang dibuat seperti ini:
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
Kolom s selalu diisi dengan nol, yang saya tidak mengerti, gambar dimasukkan dalam bentuk biner asli, format (jpg, png, gif) hilang, tetapi OsmAnd menentukannya dari isinya. Ubin dalam format berbeda dapat dikemas dalam satu database. Alih-alih $ minzoom dan $ maxzoom, itu perlu untuk mengganti batas skala ubin yang dimasukkan di basis. File database yang telah selesai, misalnya, Topo.sqlitedb, ditransfer ke smartphone atau tablet di direktori osmand / tiles. Mulai ulang OsmAnd, di Menu -> "Configure Map" -> "Top Layer" opsi baru "Topo" akan muncul - ini adalah peta dengan ubin baru kami.