في الجزء السابق ، توقفنا عن حقيقة أنه من الضروري الحصول على لون البكسل من الخريطة النقطية الأصلية بواسطة الإحداثيات الجيوديسية (خطوط
الطول والعرض) ، والتي تمت ترجمتها بالفعل إلى SC لهذه الخريطة.
يتم تعيين الإحداثيات الجيوديسية على سطح الشكل الإهليلجي وإحداثيات البكسل على المستوى ، أي أنت بحاجة إلى طريقة للانتقال من سطح محدب إلى سطح مستو. تكمن المشكلة في أن السطح المحدب (دعنا نتخيل أن نوعًا من الرسم أو شبكة إحداثيات مطبقة عليه) لا يمكن نقله إلى سطح مستوٍ دون أي تشويه. قد تكون مشوهة: الشكل (الزوايا) والمساحة والأبعاد الخطية. هناك ، بالطبع ، إمكانية ، وليست الوحيدة ، للانتقال إلى سطح مستو دون تشويه شيء واحد فقط ، لكن البقية ستشوه حتما.
من الواضح ، على المقاييس الأصغر (الكوكب بأكمله ، القارة) ، أن التشوهات تكون أكثر وضوحًا من التشوهات الأكبر (المنطقة ، المدينة) ، وعلى أكبر (مخطط منطقة صغيرة) ، لا نتحدث عنها على الإطلاق ، لأن لم يعد من الممكن تمييز سطح الشكل الإهليلجي على مثل هذه المقاييس عن المستوى المستوي.
هنا نأتي إلى مفهوم إسقاط الخريطة. لن أعطي أمثلة بصور إسقاط كرة على أسطوانة أو مخروط مع التطور اللاحق ، بالنسبة لنا من المهم الآن التعميم.
إسقاط الخريطة هو طريقة محددة رياضيًا لرسم خرائط لسطح مجسم إهليلجي على مستوى. ببساطة ، هذه صيغ رياضية لتحويل الإحداثيات الجيوديسية (خطوط الطول والعرض) إلى إحداثيات ديكارتية مسطحة - فقط ما نحتاجه.
تم اختراع مجموعة كبيرة ومتنوعة من الإسقاطات الخرائطية ، وجميعهم يجدون تطبيقًا لأغراضهم الخاصة: الحجم المتساوي (حيث يتم الحفاظ على المنطقة) للخرائط السياسية ، وخرائط التربة ، والمطابقة (التي يتم فيها الحفاظ على الشكل) - للملاحة ، متساوية البعد (الحفاظ على المقياس في الاتجاه المختار) - للكمبيوتر معالجة وتخزين صفائف البيانات الجغرافية. هناك أيضًا توقعات تجمع بين الميزات المذكورة أعلاه بنسب معينة ، وهناك توقعات مقصورة على فئة معينة تمامًا. يمكن العثور على أمثلة من الإسقاطات على ويكيبيديا في صفحة قائمة إسقاطات الخريطة.
لكل إسقاط ، يتم اشتقاق الصيغ الدقيقة ، أو في شكل مجموع سلسلة متقاربة لا نهائية ، وفي الحالة الأخيرة ، قد يكون هناك عدة خيارات. تقوم صيغ الإسقاط بتحويل خطوط الطول والعرض إلى إحداثيات ديكارتية ، وعادة ما يتم استخدام العداد كوحدة قياس في مثل هذه الإحداثيات. يمكن أحيانًا رسم شبكة مستطيلة يبلغ طولها مترًا واحدًا على خريطة (في شكل شبكة كيلومتر) ، ولكن في معظم الحالات لا يتم رسمها. لكننا نعلم الآن أنه لا يزال موجودًا بشكل ضمني. يتم حساب مقياس الخريطة ، المشار إليه على الخريطة ، فقط بالنسبة إلى حجم هذه الشبكة. يجب أن يكون مفهوماً بوضوح أن مترًا واحدًا على شبكة الإحداثيات هذه يتوافق مع متر واحد على الأرض ليس على الخريطة بأكملها ، ولكن فقط في بعض النقاط ، على طول خط معين ، أو على طول الخطوط في اتجاه معين ،في بقية النقاط أو الاتجاهات ، تظهر تشوهات و 1 متر على الأرض لا يتوافق مع متر واحد من شبكة الإحداثيات.
استطرادية صغيرة. الوظيفة من الجزء الأول من المقالة WGS84_XYZ هي على وجه التحديد تحويل الإحداثيات من WGS84 SC إلى إحداثيات مستطيلة ، ولكن ليس بالأمتار ، ولكن بالبكسل. كما ترون ، فإن الصيغة هناك بسيطة للغاية ، إذا لم يتم استخدام إسقاط Mercator على كرة ، ولكن على شكل بيضاوي ، فستكون الصيغة أكثر تعقيدًا. لهذا السبب ، من أجل تسهيل الحياة على المتصفحات ، تم اختيار مجال ، وبالتالي ترسخ إسقاط WebMercator مع مجاله ، والذي غالبًا ما يتم انتقاده بسببه.
بالنسبة لاحتياجاتي ، استخدمت مستندًا يسمى "إسقاطات الخرائط التي يستخدمها المسح الجيولوجي الأمريكي" بتنسيق pdf ، والتي يمكن العثور عليها على الإنترنت. يوفر المستند إرشادات مفصلة لكل إسقاط لتسهيل كتابة دالة تحويل بلغة برمجة.
دعنا نواصل كتابة الخوارزمية الخاصة بنا. دعنا ننفذ أحد الإسقاطات الشائعة المسماة Transverse Mercator وأحد متغيراتها يسمى إسقاط 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,
};
تتمثل إحدى ميزات إسقاط Mercator المستعرض في أنه مطابق ، أي أنه يتم الحفاظ على شكل الكائنات على الخريطة والزوايا ، بالإضافة إلى حقيقة أن المقياس محفوظ على طول خط طول مركزي واحد. يُعد الإسقاط مناسبًا للكرة الأرضية بأكملها ، ولكن التشوهات تزداد مع المسافة من خط الطول المركزي ، لذلك يتم قطع سطح الأرض بالكامل إلى شرائح ضيقة على طول خطوط الطول ، والتي تسمى مناطق ، يتم استخدام خط الطول المركزي لكل منها. أمثلة على تنفيذ مثل هذه الإسقاطات هي إسقاط Gauss-Kruger و 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) إحداثيات ديكارتية مستطيلة بالأمتار.
نحن بالفعل قريبون جدًا من إكمال الخوارزمية الخاصة بنا. علينا فقط إيجاد إحداثيات البكسل (س ، ص) في ملف الخريطة. لان إحداثيات البكسل هي أيضًا ديكارتية ، ثم يمكننا العثور عليها عن طريق التحويل التقريبي (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:
في الخوارزمية ، ظلت وظيفة العثور على الصورة الأصلية لخريطة srcimg من الإحداثيات الجيوديسية lat ، lon المحددة في SC للخريطة غير موصوفة. أعتقد أنه لن تكون هناك مشاكل معها وعدد منطقة المنطقة> 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;
}
}
عند الإدخال ، يجب تقديم ثلاث نقاط ربط على الأقل ، عند الإخراج نحصل على معلمات جاهزة. نقاط الارتساء هي نقاط معروفة لكل من إحداثيات البكسل (س ، ص) وإحداثيات إزاحة الشرق والشمال (هـ ، ن). يمكن استخدام نقاط تقاطع شبكة الكيلومترات على الخريطة الأصلية كنقاط كهذه. ماذا لو لم تكن هناك شبكة كيلومترات على الخريطة؟ ثم يمكنك تعيين الأزواج (x ، y) و (lon ، lat) ، حيث تأخذ هذه النقاط نقاط تقاطع المتوازيات وخطوط الطول ، فهي دائمًا على الخريطة. يبقى فقط التحويل (lon ، lat) إلى (e ، n) ، ويتم ذلك من خلال وظيفة التحويل للإسقاط ، وفي حالتنا ، يتم ترجمته ().
كما ترى ، هناك حاجة إلى نقاط الربط لإخبار الخوارزمية بأي جزء من سطح الأرض يصفه الملف مع صورة الخريطة. نظرًا لأن كلا نظامي الإحداثيات كانا ديكارتيين ، بغض النظر عن عدد نقاط الربط التي حددناها وبغض النظر عن بعدهما عن بعضهما البعض ، فإن التناقض على مستوى الخريطة بأكمله سيكون ضمن الخطأ في تحديد نقاط الربط. معظم الأخطاء هي أن الإسقاط أو الإسناد أو الشكل الإهليلجي الخاطئ (مع عدم تحديد المعلمات بدقة) يتم استخدامه ، ونتيجة لذلك ، فإن الإحداثيات (e ، n) عند الإخراج ليست في نظام الإحداثيات الديكارتية ، ولكن في منحني قليلاً بالنسبة إلى الديكارتي. بصريا ، يمكن تصور هذا على أنه "ورقة مجعدة". بطبيعة الحال ، فإن زيادة عدد نقاط الربط لا يحل هذه المشكلة. يمكن حل المشكلة عن طريق ضبط معلمات الإسقاط والمرجع والقطع الناقص ،في هذه الحالة ، سيسمح لك عدد كبير من نقاط الربط بتنعيم "الورقة" بشكل أكثر دقة وعدم تفويت المناطق غير الملساء.
وفي نهاية المقال سأخبرك بكيفية استخدام البلاط الجاهز في OpenLayers وبأي شكل لإعدادهم لبرنامج OsmAnd.
بالنسبة إلى OpenLayers ، تحتاج المربعات ببساطة إلى نشرها على الويب وتسميتها بحيث يمكن تمييز (z، y، x) في اسم الملف أو الدليل ، على سبيل المثال:
/tiles/topo/12_1377_2391.jpg
أو أفضل:
/ بلاط / 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
يتم تعبئة الأعمدة دائمًا بالصفر ، والتي لم أفهمها ، يتم إدخال الصورة بالشكل الثنائي الأصلي ، ويتم فقد التنسيق (jpg ، png ، gif) ، ولكن OsmAnd يحددها بمحتواها. يمكن حزم المربعات بتنسيقات مختلفة في قاعدة بيانات واحدة. بدلاً من minzoom $ و maxzoom $ ، من الضروري استبدال حدود المقياس للبلاطات التي تم إدخالها في القاعدة. يتم نقل ملف قاعدة البيانات المكتمل ، على سبيل المثال ، Topo.sqlitedb ، إلى هاتف ذكي أو جهاز لوحي في دليل osmand / البلاط. أعد تشغيل OsmAnd ، في القائمة -> "تكوين الخريطة" -> "الطبقة العليا" ، سيظهر خيار جديد "Topo" - هذه خريطة مع مربعاتنا الجديدة.