إنشاء مربعات من الخرائط النقطية

بطريقة ما شعرت بالحيرة من مسألة إنشاء خرائط مناسبة للاستخدام في OsmAnd و OpenLayers. في ذلك الوقت لم يكن لدي أي فكرة عن نظم المعلومات الجغرافية على الإطلاق ، لذلك تعاملت مع كل شيء من البداية.



سأخبرك في المقالة عن نتائج "بحثي" ، قم بتكوين خوارزمية لتحويل خريطة نقطية عشوائية إلى مربعات مفهومة للتطبيقات ، وعلى طول الطريق ، تعرف على مفاهيم مثل الشكل الإهليلجي ، والمرجع ، ونظام الإحداثيات ، والإسقاط.



أرضنا ليست على شكل كرة ، ولا حتى على شكل إهليلجي ؛ هذا الشكل المعقد يسمى الجيود. الحقيقة هي أن الكتل داخل الأرض ليست موزعة بالتساوي ، لذلك في بعض الأماكن تكون الأرض مقعرة قليلاً ، وفي أماكن أخرى تكون محدبة قليلاً. إذا أخذنا أراضي دولة أو قارة منفصلة ، فيمكن محاذاة سطحها بدقة كافية مع سطح بعض القطع الناقص ، إذا تم إزاحة مركز هذا الشكل الإهليلجي قليلاً على طول ثلاثة محاور إحداثية بالنسبة إلى مركز كتلة الأرض. يسمى هذا الشكل الإهليلجي مرجعيًا إهليلجيًا ، وهو مناسب لوصف فقط المنطقة المحلية من الأرض التي تم إنشاؤها من أجلها. على مسافات كبيرة من هذا الموقع ، يمكن أن يكون لها تناقض كبير جدًا مع سطح الأرض. يسمى الشكل الإهليلجي الذي يتزامن مركزه مع مركز كتلة الأرض بالمجسم الإهليلجي الأرضي المشترك. واضح،أن الشكل الإهليلجي المرجعي يصف منطقته المحلية من الأرض بشكل أفضل من وصف الأرض العامة ، لكن الأرض العامة مناسبة لسطح الأرض بأكمله.



لوصف الشكل الإهليلجي ، لا تكفي سوى قيمتين مستقلتين: نصف القطر الاستوائي (يُشار إليه عادةً بعلامة أ) ونصف القطر القطبي (ب) ، ولكن بدلاً من القيمة المستقلة الثانية ، يُستخدم الانكماش القطبي f = (أب) / أ عادةً. هذا هو أول شيء نحتاجه في الخوارزمية لدينا ككائن - شكل بيضاوي. بالنسبة لأجزاء مختلفة من الأرض في سنوات مختلفة ، قام باحثون مختلفون بحساب العديد من الأشكال الإهليلجية المرجعية ، ويتم تقديم المعلومات عنها في شكل بيانات: أ (بالأمتار) و 1 / و (بلا أبعاد). من الغريب أنه بالنسبة إلى الأشكال الإهليلجية الأرضية الشائعة ، هناك أيضًا العديد من المتغيرات المختلفة (مختلفة أ ، و) ، لكن الفرق ليس قويًا جدًا ، ويرجع ذلك أساسًا إلى الاختلاف في طرق تحديد 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,
};

يُظهر المثال اثنين من الأشكال الإهليلجية: WGS84 الأرضي المشترك ، المستخدم في الملاحة عبر الأقمار الصناعية ، والإهليلجي المرجعي Krasovsky ، المطبق على أراضي أوروبا وآسيا.



ضع في اعتبارك نقطة أخرى مثيرة للاهتمام ، والحقيقة هي أن شكل الأرض بطيء ، لكنه يتغير ، لذا فإن الشكل الإهليلجي ، الذي يصف السطح اليوم جيدًا ، قد يكون بعيدًا عن الواقع خلال مائة عام. هذا لا علاقة له بالمجسم الإهليلجي الأرضي المشترك ، منذ ذلك الحين الانحرافات ضمن نفس الخطأ ، ولكنها ذات صلة بالمجسم الإهليلجي المرجعي. هنا نأتي إلى مفهوم آخر - مرجع. الإسناد عبارة عن مجموعة من معلمات الشكل الإهليلجي (أ ، و) ، إزاحتها داخل الأرض (للإشارة الإهليلجية) ، مثبتة في لحظة معينة من الزمن. بتعبير أدق ، قد لا يصف المرجع بالضرورة الشكل الإهليلجي ، بل يمكن أن يكون معلمات لشخص أكثر تعقيدًا ، على سبيل المثال ، شبه منحرف.



من المؤكد أن السؤال قد نشأ بالفعل: كيف تنتقل من شكل بيضاوي أو مرجع إلى آخر؟ للقيام بذلك ، يجب أن يكون لكل شكل بيضاوي نظام إحداثيات جيوديسي: خطوط الطول والعرض (فاي ، لامدا) ، ويتم الانتقال عن طريق ترجمة الإحداثيات من نظام إحداثيات إلى آخر.

هناك العديد من الصيغ لتحويل الإحداثيات. يمكنك أولاً ترجمة الإحداثيات الجيوديسية في نظام إحداثي واحد إلى إحداثيات ثلاثية الأبعاد X و Y و Z ، وإجراء عمليات النقل والدوران معها ، ثم تحويل الإحداثيات ثلاثية الأبعاد الناتجة إلى إحداثيات جيوديسية في نظام إحداثيات آخر. يمكنك القيام بذلك مباشرة. لان جميع الصيغ عبارة عن سلسلة متقاربة لانهائية ، ثم تقتصر عادةً على عدد قليل من أعضاء السلسلة لتحقيق الدقة المطلوبة. كمثال ، سوف نستخدم تحويلات هيلميرت ، هذه تحويلات باستخدام الانتقال إلى إحداثيات ثلاثية الأبعاد ، وهي تتكون من المراحل الثلاث الموضحة أعلاه. للتحولات ، بالإضافة إلى اثنين من الأشكال البيضاوية ، تحتاج إلى 7 معلمات أخرى: ثلاث تحولات على ثلاثة محاور ، وثلاث زوايا دوران ، وعامل مقياس. كما قد تتخيل ، يمكن استخلاص جميع المعلمات من البيانات.لكن في الخوارزمية لن نستخدم مثل هذا الكائن كمرجع ، ولكن بدلاً من ذلك سنقدم كائنًا انتقاليًا من نظام إحداثي إلى آخر ، والذي سيحتوي على: إشارات إلى اثنين من الأشكال البيضاوية و 7 معلمات تحويل. سيكون هذا هو العنصر الثاني في خوارزمية لدينا.

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},
};

يوضح المثال معلمات التحويل من نظام إحداثيات Pulkovo 1942 إلى نظام إحداثيات 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 وفي الويب ، على سبيل المثال ، في خرائط Google وفي OpenStreetMap. WebMercator هو إسقاط مركاتور مبني على كرة. الإحداثيات في هذا الإسقاط هي إحداثيات البكسل X ، Y وتعتمد على مقياس Z ، بالنسبة لمقياس صفري ، يتم وضع سطح الأرض بالكامل (حتى حوالي 85 درجة من خط العرض) على قطعة واحدة 256 × 256 بكسل ، تتغير إحداثيات 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);
}

وفي نهاية الجزء الأول من المقالة ، يمكننا بالفعل رسم بداية الخوارزمية الخاصة بنا لإنشاء مربع. تتم معالجة كل بلاطة بحجم 256 × 256 بكسل بثلاث قيم: x ، y ، z ، العلاقة مع الإحداثيات X ، Y ، Z بسيطة جدًا: x = (int) (X / 256) ؛ y = (int) (Y / 256) ؛ ض = ض ؛

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. ستكون هذه هي المقالة التالية ، التي سنتحدث فيها عن الإسقاطات الجيوديسية والتحولات الأفينية.