كيف تعمل الخرائط الرقمية؟ الإحداثيات والمربعات Tiles

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

هذا المقال هو أساس مرحلة الخرائط، وما بعده من استضافة الخرائط والبحث عن العناوين وحساب المسارات يفترض أنك تعرف ما فيه. وسوف نناقش في هذا المقال ما يلي:

• خط العرض Latitude وخط الطول Longitude، ونظام الإحداثيات WGS84 ورقمه EPSG:4326.
• الإسقاط Projection، وإسقاط Web Mercator ورقمه EPSG:3857 الذي تستخدمه أغلب خرائط الويب، ولماذا يكبر المساحات قرب القطبين.
• المربعات Tiles ومستويات التقريب Zoom Levels، وحساب عنوان المربع z/x/y بالأرقام لمدينة الرياض.
• طبقات بيانات الخريطة: الأرض والماء والطرق والمباني والأماكن والنصوص، والأسماء العربية.
• دقة الإحداثيات: كم متراً يمثل كل رقم بعد الفاصلة، ثم GeoHash وH3 باختصار.
• تجربة عملية: نحسب المربعات بلغة Python، ونفتح مربعاً حقيقياً من ملف PMTiles ونقرأ ما فيه، ونرسم حدود المربعات على خريطة الرياض في المتصفح، ونقارن حساب المسافة بثلاث طرق.

خط العرض وخط الطول

أي نقطة على الأرض تحدد برقمين، وهما زاويتان وليسا مسافتين:

• خط العرض Latitude: الزاوية شمال خط الاستواء أو جنوبه، من 90 جنوباً إلى 90 شمالاً، فخط الاستواء 0، والرياض قرابة 24.64 شمالاً.
• خط الطول Longitude: الزاوية شرق خط جرينتش Prime Meridian أو غربه، من 180 غرباً إلى 180 شرقاً، والرياض قرابة 46.72 شرقاً.

وتكتب الأرقام في البرمجة بالدرجات العشرية Decimal Degrees، والجنوب والغرب بإشارة سالبة، فنقطة الرياض في OpenStreetMap هي 24.6389160, 46.7160104. وانتبه للترتيب، فهو أكثر خطأ ستقع فيه: الإنسان يقول «خط العرض ثم الطول»، وأغلب المعايير والمكتبات البرمجية، مثل GeoJSON وMapLibre وPostGIS، تكتب الطول أولاً [46.716, 24.639]، لأنه يقابل المحور الأفقي x. فإذا ظهرت نقاطك في المحيط الهندي بدل الرياض، فأول ما تفحصه هو الترتيب.

ولكن الزاوية تحتاج شكلاً تقاس عليه، والأرض ليست كرة تامة وإنما مفلطحة قليلاً عند القطبين، لذلك تحتاج نموذجاً رياضياً للشكل اسمه المرجع الجيوديسي Geodetic Datum، وأشهرها WGS84 الذي يستخدمه GPS، ورقمه في سجل EPSG هو 4326. فعندما ترى EPSG:4326 في أي أداة فالمعنى: خط عرض وخط طول بالدرجات على نموذج WGS84.

الإسقاط: كيف تصبح الكرة مستطيلاً؟

الشاشة مسطحة والأرض كرة، ولا توجد طريقة لفرد قشرة برتقالة على طاولة دون أن تتمزق أو تتمدد، وهذه هي مشكلة الإسقاط Projection: دالة رياضية تحول خط العرض والطول إلى x وy على سطح مستو، وكل إسقاط يضحي بشيء ليحفظ شيئاً آخر: المساحة أو الأشكال أو المسافات أو الاتجاهات، ولا يحفظها كلها.

والإسقاط الذي تستخدمه أغلب خرائط الويب هو Web Mercator، ورقمه EPSG:3857، واسمه الرسمي في السجل WGS 84 / Pseudo-Mercator. وسبب اختياره عملي: يحفظ الأشكال المحلية والزوايا، فالتقاطع القائم في الشارع يبقى قائماً على الشاشة، والشمال دائماً للأعلى، والعالم كله يصبح مربعاً واحداً يسهل تقسيمه إلى مربعات أصغر. والثمن أن المقياس يكبر كلما ابتعدت عن خط الاستواء، فالطول عند خط عرض φ يكبر بالمعامل 1 / cos(φ) والمساحة بمربعه، ولهذا يتوقف الإسقاط عند خط عرض 85.0511 تقريباً، لأن المعامل عند القطب يصبح لا نهائياً، وهذا الحد مذكور في توثيق OpenStreetMap للمربعات.

وهذا يفسر جرينلاند: هي قرب خط عرض 72، فكل كيلومتر مربع فيها يرسم أكبر بقرابة 10.5 مرة من كيلومتر مربع عند خط الاستواء، وإفريقيا أغلبها قرب خط الاستواء فترسم بحجم قريب من حقيقتها، فتبدوان متقاربتين على الخريطة وإفريقيا في الحقيقة أكبر بكثير. والرياض عند خط عرض 24.64 معاملها 1.10 للطول و1.21 للمساحة، وسوف نحسب هذه الأرقام في التجربة.

والقاعدة العملية التي تخرج بها: Web Mercator ممتاز للعرض والتنقل، وسيئ لحساب المساحات والمسافات، فلا تحسب مساحة أرض أو طول طريق من x وy على الخريطة، وإنما من خط العرض والطول على النموذج الجيوديسي، كما سوف نرى في التجربة.

المربعات ومستويات التقريب

بعد الإسقاط أصبح العالم مربعاً واحداً، والآن نقسمه، فقد شرحنا في خرائطك على سيرفرك أن الخريطة على الويب لا تحمل كصورة واحدة، وإنما كمربعات Tiles في هرم من مستويات التقريب، والفرق بين المربع النقطي Raster والمتجه Vector، ولن نكرر ذلك هنا. ما نضيفه هو الحساب نفسه:

• في المستوى z يقسم العالم إلى 2^z عموداً و2^z صفاً، أي 4^z مربعاً: مربع واحد في المستوى 0، وأربعة في المستوى 1، و1,024 في المستوى 5، وأكثر من مليار في المستوى 15.
• العمود x يبدأ من 0 عند خط الطول 180 غرباً ويزيد شرقاً، والصف y يبدأ من 0 في الأعلى عند خط عرض 85.0511 شمالاً ويزيد جنوباً، أي عكس اتجاه خط العرض.
• عنوان المربع هو z/x/y، ومنه يأتي رابط المربع في أي سيرفر خرائط، مثل /tiles/12/2579/1758.png.

والمعادلة من توثيق OpenStreetMap:

n     = 2 ^ zoom
xtile = n × (lon + 180) / 360
ytile = n × (1 − ln(tan(lat) + sec(lat)) / π) / 2

والجزء الصحيح من الناتج هو رقم المربع، والكسر هو مكان النقطة داخله، فإذا كان x = 20.15 فالنقطة في العمود 20 وعلى بعد 15% من حافته اليسرى. وسوف نرى ذلك بالصورة في التجربة. والصورة التالية تبين الترقيم في المستوى 2 ومكان الرياض في المستوى 5:

وكم متراً يغطي البكسل الواحد؟ عرض الأرض عند خط الاستواء قرابة 40,075 كيلومتراً، يقسم على عدد البكسلات في المستوى، وهو 256 مضروباً في 2^z للمربعات التقليدية بعرض 256 بكسل، ثم يضرب في cos(lat) بسبب الإسقاط. ففي الرياض يغطي البكسل قرابة 139 متراً في المستوى 10، و4.3 متر في المستوى 15، وهذا ما يحدد في أي مستوى يبدأ ظهور المباني أو أسماء الشوارع.

طبقات الخريطة: ما الذي يرسم فوق ماذا؟

الخريطة التي تراها طبقات فوق بعضها، والمربع المتجه يحمل كل طبقة باسمها، والمتصفح يرسمها بالترتيب الذي يحدده ملف التنسيق Style:

• الأساس Base: الأرض earth والماء water، ثم استخدامات الأرض landuse مثل الحدائق والمناطق الصناعية.
• الشبكات Networks: الطرق roads بأنواعها، والسكك الحديدية.
• المباني Buildings: مضلعات تظهر في المستويات القريبة فقط.
• الأماكن Places ونقاط الاهتمام Points of Interest واختصاراً POI: أسماء المدن والأحياء، والمستشفيات والمطاعم والمساجد، وكل نقطة لها نوع kind وأسماء بعدة لغات.
• النصوص Labels: ترسم أخيراً فوق كل شيء، والمكتبة تقرر أي نص يظهر وأيها يخفى حتى لا تتزاحم.

ولأن الأسماء بيانات داخل المربع المتجه وليست بكسلات مرسومة، فنفس المربع يعرض العربية أو الإنجليزية بحسب ملف التنسيق، والاسم العربي يأتي من الوسم name:ar في OpenStreetMap. ورسم النص العربي متصلاً ومن اليمين إلى اليسار كان مشكلة قديمة في مكتبات الخرائط، وقد اختبرناه في النصوص العربية على الخريطة، والخلاصة أن MapLibre من الإصدار 6.9.0 ترسمه صحيحاً دون إضافة، وأن الاسم الناقص بالعربية يصلح في البيانات لا في الكود.

دقة الإحداثيات: كم رقماً بعد الفاصلة تحتاج؟

درجة خط العرض قرابة 111 كيلومتراً في أي مكان، أما درجة خط الطول فتصغر كلما ابتعدت عن خط الاستواء بنفس معامل cos(lat)، ففي الرياض قرابة 101 كيلومتر. ومن هنا تعرف ماذا يعني كل رقم بعد الفاصلة:

والقاعدة العملية: ستة أرقام بعد الفاصلة تكفي أي تطبيق تقريباً، وخمسة تكفي توصيل الطلبات، وكل رقم بعدها يخزن ضجيجاً لا دقة. والعكس مهم في الخصوصية: إذا عرضت مواقع العملاء أو المندوبين في تقرير، فقرب الأرقام إلى رقمين أو ثلاثة حتى لا تكشف المنزل بعينه. ولا تخزن الإحداثيات كنوع float بدقة مفردة، لأن دقته قرابة 7 أرقام معنوية فقط، فخط طول مثل 46.7160104 يخزن بخطوات قرابة 0.0000038 درجة، أي قرابة 40 سنتيمتراً، فما بعد الرقم الخامس بعد الفاصلة لا يحفظ بدقة، وقد شرحنا السبب في الأنظمة العددية للمبرمجين، واستخدم double أو نوعاً جغرافياً مثل geography في PostGIS.

GeoHash وH3: الإحداثيات كخلايا

أحياناً لا تحتاج النقطة نفسها، وإنما الخلية التي هي فيها: كم طلباً في كل منطقة؟ أي السائقين قريب من هذا العميل؟ وهنا تأتي أنظمة تقسم الأرض إلى خلايا لها أسماء:

• GeoHash: يحول النقطة إلى نص قصير، كل حرف إضافي فيه يقسم الخلية إلى خلايا أصغر، فالنقاط التي تشترك في بداية النص قريبة من بعضها، وهذا يجعل البحث عن القريب بحثاً عن بادئة نص في فهرس عادي، كما في الدالة ST_GeoHash في PostGIS. وخلاياه مستطيلات تختلف أبعادها بحسب خط العرض.
• H3: نظام من Uber برخصة Apache 2 يقسم الأرض إلى خلايا سداسية Hexagons، في 16 مستوى من 0 إلى 15، وفي كل مستوى 12 خلية خماسية فقط لا يمكن تجنبها. وميزة السداسي أن كل خلية لها ستة جيران على نفس المسافة تقريباً، فيصلح لتحليل الطلب وتسعير الرحلات والخرائط الحرارية.

التجربة العملية

سوف نستخدم Python داخل Container مع GeographicLib للحساب الدقيق على نموذج WGS84، وgo-pmtiles لقص مربعات الرياض من بناء Protomaps اليومي، وصفحة MapLibre من مختبر خرائطك على سيرفرك مع تعديل صغير يرسم حدود المربعات. والنقطة في كل التجربة هي عقدة الرياض في OpenStreetMap، رقم 315358390.

حساب عنوان المربع للرياض

ملف tiles.py:

# Tile math for one point: the Riyadh node in OpenStreetMap (node 315358390)
import math
lat, lon = 24.6389160, 46.7160104

def tile(lat, lon, z):
    n = 2 ** z
    x = (lon + 180) / 360 * n
    lat_r = math.radians(lat)
    y = (1 - math.asinh(math.tan(lat_r)) / math.pi) / 2 * n
    return x, y

print("zoom  tiles_in_world   x            y            tile")
for z in (0, 1, 5, 10, 12, 15):
    x, y = tile(lat, lon, z)
    print(f"{z:>4}  {4**z:>14,}  {x:>11.4f}  {y:>11.4f}  {z}/{int(x)}/{int(y)}")

docker run --rm --name lab-mapsb-py -v "$PWD:/lab" -w /lab python:3.13.16-slim python tiles.py

والمخرج سوف يكون كما يلي:

zoom  tiles_in_world   x            y            tile
   0               1       0.6298       0.4293  0/0/0
   1               4       1.2595       0.8587  1/1/0
   5           1,024      20.1525      13.7391  5/20/13
  10       1,048,576     644.8811     439.6504  10/644/439
  12      16,777,216    2579.5244    1758.6014  12/2579/1758
  15   1,073,741,824   20636.1951   14068.8113  15/20636/14068

في المخرج أعلاه لاحظ التالي:

• math.asinh(math.tan(lat)) هي نفس ln(tan(lat) + sec(lat)) في المعادلة، بصيغة أقصر.
• في المستوى 1 العالم أربعة مربعات، والرياض في المربع 1/1/0: العمود 1 لأنها شرق جرينتش، والصف 0 لأنها شمال خط الاستواء.
• كل مستوى يضاعف x وy تقريباً: 20.15 في المستوى 5 تصبح 644.88 في المستوى 10، أي ضرب في 2 خمس مرات، وهذا ما يجعل المربع في مستوى ما يحتوي أربعة مربعات من المستوى الذي بعده.
• عدد المربعات في المستوى 15 أكثر من مليار، لذلك لا أحد يرسم المربعات النقطية للعالم كله مسبقاً، وإنما يقص المنطقة التي يحتاجها أو يستخدم المربعات المتجهة.

مربع حقيقي: ماذا يوجد داخل 12/2579/1758؟

نقص مستطيل الرياض من بناء Protomaps بتاريخ 20261008 حتى المستوى 14، كما فعلنا في مختبر خرائطك على سيرفرك، ثم نستخرج المربع الذي حسبناه للمستوى 12:

mkdir -p web/tiles
docker run --rm --name lab-mapsb-pmtiles -v "$PWD/web/tiles:/out" protomaps/go-pmtiles:v1.31.2 \
  extract https://build.protomaps.com/20261008.pmtiles /out/riyadh.pmtiles --bbox=46.45,24.45,47.05,25.00 --maxzoom=14
docker run --rm --name lab-mapsb-pmtiles -v "$PWD/web/tiles:/out" protomaps/go-pmtiles:v1.31.2 \
  tile /out/riyadh.pmtiles 12 2579 1758 > tile-12-2579-1758.mvt

والمخرج من أمر القص سوف يكون كما يلي (اختصرنا شريط التقدم):

extract.go:450: fetching 1108 tiles, 57 chunks, 27 requests
extract.go:606: Completed in 32.584997383s with 4 download threads (34.00337841296619 tiles/s).
extract.go:611: Extract required 37 total requests.
extract.go:612: Extract transferred 9.4 MB (overfetch 0.05) for an archive size of 9.1 MB

ثم نفتح المربع بمكتبة mapbox-vector-tile ونعد ما في كل طبقة. ملف layers.py:

# Open one vector tile (z12/2579/1758) and list what is inside it
import gzip, mapbox_vector_tile
raw = gzip.decompress(open("tile-12-2579-1758.mvt", "rb").read())
tile = mapbox_vector_tile.decode(raw)
print(f"compressed 99121 bytes, uncompressed {len(raw)} bytes")
for name, layer in tile.items():
    kinds = sorted({f["geometry"]["type"] for f in layer["features"]})
    print(f"{name:<10} {len(layer['features']):>5} features  extent {layer['extent']}  {', '.join(kinds)}")
pois = [f["properties"] for f in tile["pois"]["features"] if "name:ar" in f["properties"]]
print("\nfirst POIs with an Arabic name:")
for p in pois[:4]:
    print({k: p[k] for k in ("kind", "name", "name:ar", "name:en") if k in p})

docker run --rm --name lab-mapsb-py -v "$PWD:/lab" -w /lab python:3.13.16-slim sh -c 'pip install -q --root-user-action=ignore mapbox-vector-tile==2.2.0 2>/dev/null; python layers.py'

والمخرج سوف يكون كما يلي:

compressed 99121 bytes, uncompressed 131382 bytes
buildings      2 features  extent 4096  MultiPolygon
earth          3 features  extent 4096  MultiLineString, Polygon
landuse      318 features  extent 4096  MultiPolygon, Polygon
places        49 features  extent 4096  Point
pois          57 features  extent 4096  Point
roads        155 features  extent 4096  LineString, MultiLineString
water          4 features  extent 4096  MultiLineString, MultiPolygon, Polygon

first POIs with an Arabic name:
{'kind': 'hospital', 'name': 'مستشفى الملك فيصل التخصصي', 'name:ar': 'مستشفى الملك فيصل التخصصي', 'name:en': 'King Faisal Specialist Hospital'}
{'kind': 'castle', 'name': 'قصر الناصرية', 'name:ar': 'قصر الناصرية', 'name:en': 'Qasr an Nasriyah'}
{'kind': 'administrative', 'name': 'الفاخرية', 'name:ar': 'الفاخرية', 'name:en': 'Al Fakhiriyah District'}
{'kind': 'administrative', 'name': 'البديعة', 'name:ar': 'البديعة', 'name:en': 'Al Badiah District'}

في المخرج أعلاه لاحظ التالي:

• المربع كله 99 كيلوبايت مضغوطاً و128 كيلوبايت دون ضغط، وفيه كل ما تحتاجه لرسم وسط الرياض في المستوى 12 بأي لغة وأي ألوان.
• الطبقات السبع هي نفسها طبقات الخريطة التي شرحناها، وكل عنصر فيها هندسة Geometry ونوع، فالطرق خطوط والأماكن نقاط والمباني مضلعات.
• extent 4096 تعني أن الإحداثيات داخل المربع أعداد صحيحة من 0 إلى 4096، وليست خطوط عرض وطول، فالمكتبة تحولها إلى بكسلات بحسب حجم المربع على الشاشة، وهذا من أسباب صغر حجمه، كما في مواصفة Mapbox Vector Tile.
• المباني عنصران فقط في المستوى 12، لأن البناء لا يضع المباني إلا في المستويات القريبة، وكل نقطة اهتمام فيها name وname:ar وname:en، فالتنسيق يختار منها.

نرى حدود المربعات على الخريطة

نضيف إلى مجلد web ملفات MapLibre وpmtiles والخطوط كما في ملفات الموقع في المقال السابق، وبنفس إعداد nginx، ثم صفحة جديدة web/tiles.html تفعل الخيار showTileBoundaries الذي يرسم حدود كل مربع وعنوانه وحجمه، وتسجل عنوان كل مربع تحمله الخريطة:

<!doctype html>
<html lang="ar" dir="rtl">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>مربعات الرياض</title>
<link rel="stylesheet" href="/lib/maplibre-gl.css">
<style>html, body, #map { margin: 0; height: 100%; }</style>
<script type="importmap">{ "imports": { "fflate": "/lib/fflate.mjs" } }</script>
</head>
<body>
<div id="map"></div>
<script type="module">
import * as maplibregl from '/lib/maplibre-gl.mjs';
import { Protocol } from '/lib/pmtiles.mjs';
import { layers, namedFlavor } from '/lib/basemaps.mjs';

maplibregl.addProtocol('pmtiles', new Protocol().tile);
const q = new URLSearchParams(location.search);
const riyadh = [46.7160104, 24.6389160];          // the Riyadh node in OpenStreetMap
window.map = new maplibregl.Map({
  container: 'map',
  center: riyadh,
  zoom: Number(q.get('z') || 12),
  style: {
    version: 8,
    glyphs: '/assets/fonts/{fontstack}/{range}.pbf',
    sprite: `${location.origin}/assets/sprites/v4/light`,
    sources: { protomaps: { type: 'vector', url: 'pmtiles:///tiles/riyadh.pmtiles',
      attribution: '<a href="https://www.openstreetmap.org/copyright">© OpenStreetMap</a>' } },
    layers: layers('protomaps', namedFlavor('light'), { lang: 'ar' })
  }
});
map.showTileBoundaries = true;                    // draw each tile's edge and its z/x/y
new maplibregl.Marker({ color: '#E95420' }).setLngLat(riyadh).addTo(map);
// log every tile the map loads, as z/x/y
window.tiles = new Set();
map.on('sourcedata', (e) => {
  if (e.coord) { const c = e.coord.canonical; window.tiles.add(`${c.z}/${c.x}/${c.y}`); }
});
</script>
</body>
</html>

وملف compose.yaml بمنفذ مختلف عن المقال السابق:

name: lab-mapsb
services:
  web:
    image: nginx:1.29.8-alpine
    container_name: lab-mapsb-web
    volumes:
      - ./web:/usr/share/nginx/html:ro
      - ./nginx.conf:/etc/nginx/conf.d/default.conf:ro
    ports:
      - "127.0.0.1:18761:80"

ثم نشغله ونفتح http://127.0.0.1:18761/tiles.html?z=5، ثم z=10 وz=12 وz=14، ونقرأ قائمة المربعات من Console بالأمر JSON.stringify([...window.tiles].sort()):

{"zoom":5,"tiles":["5/20/13"]}
{"zoom":10,"tiles":["10/644/438","10/644/439","10/644/440","10/645/438","10/645/439","10/645/440"]}
{"zoom":12,"tiles":["12/2578/1757","12/2578/1758","12/2578/1759","12/2579/1757","12/2579/1758","12/2579/1759","12/2580/1757","12/2580/1758","12/2580/1759"]}
{"zoom":14,"tiles":["14/10316/7033","14/10316/7034","14/10316/7035","14/10317/7033","14/10317/7034","14/10317/7035","14/10318/7033","14/10318/7034","14/10318/7035","14/10319/7033","14/10319/7034","14/10319/7035"]}

المربع الذي حسبناه في Python لكل مستوى موجود في القائمة: 5/20/13 و10/644/439 و12/2579/1758، والمتصفح طلب حوله المربعات التي تغطي الشاشة فقط. والصورة التالية للمستوى 5، والملف الذي قصصناه فيه مربع واحد فقط في هذا المستوى، فظهر وحده وحوله رمادي:

ولاحظ مكان العلامة: حسبنا x = 20.1525 وy = 13.7391، أي 15% من عرض المربع من اليسار و74% من ارتفاعه من الأعلى، وهذا بالضبط مكانها في الصورة. والصورة التالية للمستوى 12، والمربع 12/2579/1758 في الوسط بحجم 128 كيلوبايت، وهو نفس الرقم الذي حصلنا عليه بعد فك الضغط في Python:

ولاحظ في الصور أن المربعات على الشاشة أكبر من 256 بكسل، لأن MapLibre يرسم المربعات المتجهة بعرض 512 بكسل، فمستوى التقريب الظاهر فيه يقابل مستوى أعلى بواحد في المربعات النقطية التقليدية.

المسافة: Haversine مقابل الحساب المسطح

الآن السؤال الذي يقع فيه كثير من المطورين: كم المسافة بين الرياض وجدة؟ نحسبها بأربع طرق: على نموذج WGS84 بمكتبة GeographicLib وهي المرجع، وبمعادلة Haversine على كرة، وبالطريقة البسيطة أي فيثاغورس على الدرجات كأنها شبكة مسطحة، وبفيثاغورس على x وy في Web Mercator. والمدن الثلاث من نقاطها في OpenStreetMap. ملف geo.py:

# Distances, Mercator distortion and coordinate precision, compared with the WGS84 ellipsoid (geographiclib)
import math
from geographiclib.geodesic import Geodesic

R = 6371008.8  # mean Earth radius in metres
cities = {"Riyadh": (24.6389160, 46.7160104), "Jeddah": (21.5504432, 39.1742363), "Dammam": (26.4367824, 50.1039991)}

def haversine(a, b):
    p1, p2 = math.radians(a[0]), math.radians(b[0])
    dp, dl = p2 - p1, math.radians(b[1] - a[1])
    h = math.sin(dp / 2) ** 2 + math.cos(p1) * math.cos(p2) * math.sin(dl / 2) ** 2
    return 2 * R * math.asin(math.sqrt(h))

def degrees_pythagoras(a, b):          # the simple way: treat degrees as a flat grid
    return math.hypot(b[0] - a[0], b[1] - a[1]) * 111_195

def mercator_xy(lat, lon):             # EPSG:3857 metres
    return math.radians(lon) * 6378137, math.log(math.tan(math.pi / 4 + math.radians(lat) / 2)) * 6378137

def mercator_planar(a, b):             # straight line on the Web Mercator plane
    (x1, y1), (x2, y2) = mercator_xy(*a), mercator_xy(*b)
    return math.hypot(x2 - x1, y2 - y1)

print("route            ellipsoid_km  haversine_km  err    flat_degrees_km  err    mercator_plane_km  err")
for x, y in (("Riyadh", "Jeddah"), ("Riyadh", "Dammam")):
    a, b = cities[x], cities[y]
    g = Geodesic.WGS84.Inverse(*a, *b)["s12"] / 1000
    h, f, m = haversine(a, b) / 1000, degrees_pythagoras(a, b) / 1000, mercator_planar(a, b) / 1000
    pct = lambda v: f"{(v - g) / g * 100:+.2f}%"
    print(f"{x}->{y:<9} {g:>10.1f}  {h:>12.1f}  {pct(h):>6}  {f:>12.1f}  {pct(f):>8}  {m:>14.1f}  {pct(m):>8}")

print("\nWeb Mercator scale at latitude (length factor, area factor)")
for name, lat in (("Equator", 0), ("Riyadh", 24.64), ("Greenland centre", 72), ("Max latitude", 85.0511)):
    k = 1 / math.cos(math.radians(lat))
    print(f"{name:<17} lat {lat:>7.2f}  x{k:6.2f}  area x{k*k:7.2f}")

print("\nMetres per pixel (256 px tiles) at Riyadh's latitude")
for z in (0, 5, 10, 12, 15, 18):
    print(f"z{z:<3} {40075016.686 * math.cos(math.radians(24.64)) / (256 * 2**z):>12,.2f} m")

print("\nCoordinate precision: size of the last decimal place")
for d in range(1, 8):
    step = 10 ** -d
    lat_m = Geodesic.WGS84.Inverse(24.64, 46.71, 24.64 + step, 46.71)["s12"]
    lon_m = Geodesic.WGS84.Inverse(24.64, 46.71, 24.64, 46.71 + step)["s12"]
    print(f"{d} decimals  {step:<10g} north-south {lat_m:>12,.3f} m   east-west at Riyadh {lon_m:>12,.3f} m")

docker run --rm --name lab-mapsb-py -v "$PWD:/lab" -w /lab python:3.13.16-slim sh -c 'pip install -q --root-user-action=ignore geographiclib==2.0 2>/dev/null; python geo.py'

والمخرج سوف يكون كما يلي:

route            ellipsoid_km  haversine_km  err    flat_degrees_km  err    mercator_plane_km  err
Riyadh->Jeddah         844.8         844.2  -0.07%         906.2    +7.27%           919.0    +8.79%
Riyadh->Dammam         394.5         394.3  -0.03%         426.5    +8.12%           437.5   +10.92%

Web Mercator scale at latitude (length factor, area factor)
Equator           lat    0.00  x  1.00  area x   1.00
Riyadh            lat   24.64  x  1.10  area x   1.21
Greenland centre  lat   72.00  x  3.24  area x  10.47
Max latitude      lat   85.05  x 11.59  area x 134.37

Metres per pixel (256 px tiles) at Riyadh's latitude
z0     142,289.05 m
z5       4,446.53 m
z10        138.95 m
z12         34.74 m
z15          4.34 m
z18          0.54 m

Coordinate precision: size of the last decimal place
1 decimals  0.1        north-south   11,076.829 m   east-west at Riyadh   10,124.224 m
2 decimals  0.01       north-south    1,107.676 m   east-west at Riyadh    1,012.422 m
3 decimals  0.001      north-south      110.768 m   east-west at Riyadh      101.242 m
4 decimals  0.0001     north-south       11.077 m   east-west at Riyadh       10.124 m
5 decimals  1e-05      north-south        1.108 m   east-west at Riyadh        1.012 m
6 decimals  1e-06      north-south        0.111 m   east-west at Riyadh        0.101 m
7 decimals  1e-07      north-south        0.011 m   east-west at Riyadh        0.010 m

في المخرج أعلاه لاحظ التالي:

• Haversine أخطأ بأقل من 0.1%، أي قرابة 600 متر في 845 كيلومتراً، لأنه يفترض كرة تامة بنصف قطر متوسط، وهذا يكفي لترتيب المطاعم الأقرب أو تقدير المسافة لعميل. أما إذا كنت تحسب مسافات للفوترة أو للمساحة الرسمية فاستخدم الحساب على النموذج الجيوديسي، أو ST_Distance على نوع geography في PostGIS.
• الطريقة البسيطة على الدرجات أخطأت بأكثر من 7%، أي 61 كيلومتراً، والسبب أنها تعامل درجة الطول كأنها 111 كيلومتراً مثل درجة العرض، وفي الرياض هي قرابة 101 كيلومتر فقط.
• الخط المستقيم على Web Mercator أسوأ، لأن الإسقاط يكبر كل شيء بمعامل 1.10 في هذه المنطقة كما في الجزء الثاني من المخرج، ويزيد الخطأ كلما اتجهت شمالاً، فلا تحسب مسافة من إحداثيات EPSG:3857 أبداً.
• المسافة هنا خط مستقيم على سطح الأرض، وليست مسافة الطريق بالسيارة، وهذه نحسبها بمحركات المسارات Routing في مقال لاحق في نفس المرحلة.

وآخر خطوة: نفس نقطة الرياض كخلايا GeoHash وH3. ملف cells.py:

# The same Riyadh point as a GeoHash text and as H3 hexagon cells
import math, geohash2, h3
lat, lon = 24.6389160, 46.7160104
for n in (4, 6, 8):
    g = geohash2.encode(lat, lon, precision=n)
    _, _, dlat, dlon = geohash2.decode_exactly(g)          # half the cell size, in degrees
    h = 2 * dlat * 111_320
    w = 2 * dlon * 111_320 * math.cos(math.radians(lat))
    print(f"geohash {n} chars: {g:<9} cell about {w:,.0f} x {h:,.0f} m")
for r in (5, 7, 9):
    c = h3.latlng_to_cell(lat, lon, r)
    print(f"h3 res {r}: {c}  area {h3.cell_area(c, unit='km^2'):.4f} km2  neighbours {len(h3.grid_disk(c, 1)) - 1}")

docker run --rm --name lab-mapsb-py -v "$PWD:/lab" -w /lab python:3.13.16-slim sh -c 'pip install -q --root-user-action=ignore h3==4.3.1 geohash2==1.1 2>/dev/null; python cells.py'

geohash 4 chars: th3h      cell about 35,573 x 19,568 m
geohash 6 chars: th3hpj    cell about 1,112 x 611 m
geohash 8 chars: th3hpj3s  cell about 35 x 19 m
h3 res 5: 855355d3fffffff  area 245.8690 km2  neighbours 6
h3 res 7: 875355d06ffffff  area 5.0199 km2  neighbours 6
h3 res 9: 895355d0637ffff  area 0.1024 km2  neighbours 6

لاحظ أن كل GeoHash أطول يبدأ بالأقصر منه، th3h ثم th3hpj ثم th3hpj3s، فالبحث عن كل النقاط في الخلية th3hpj هو بحث عن النصوص التي تبدأ بهذه الأحرف. وفي H3 مساحة الخلية في المستوى 7 قرابة 5 كيلومترات مربعة، ولها ستة جيران، ومعرفها رقم بصيغة ست عشرية Hexadecimal يخزن في عمود عادي.

تنظيف التجربة

كل Container من Python وgo-pmtiles كان بالخيار --rm، فيبقى سيرفر الصفحة، فنحذفه مع شبكته من داخل مجلد المختبر:

docker compose down

ثم ال Images والملفات إذا لم تعد تحتاجها:

docker image rm protomaps/go-pmtiles:v1.31.2 python:3.13.16-slim
rm -rf web tile-12-2579-1758.mvt

أخطاء شائعة

• عكس ترتيب الإحداثيات: GeoJSON وMapLibre وPostGIS تكتب الطول أولاً، وأغلب الناس والتطبيقات تكتب العرض أولاً، فإذا ظهرت النقطة في مكان غريب فابدأ من هنا.
• حساب المسافة أو المساحة على Web Mercator أو على الدرجات مباشرة: رأينا خطأ 7% إلى 11% بين مدن المملكة، ويزيد كلما اتجهت شمالاً، فاستخدم Haversine أو الحساب الجيوديسي.
• تخزين الإحداثيات بدقة مفردة float: تضيع الأمتار الأخيرة، فاستخدم double أو نوعاً جغرافياً.
• تخزين أو عرض أرقام أكثر مما تحتاج: عشرة أرقام بعد الفاصلة لا تعني دقة، وعرض ستة أرقام لموقع منزل عميل في تقرير يكشفه بالمتر.
• الخلط بين EPSG:4326 وEPSG:3857: الأول درجات والثاني أمتار على المسقط، فإذا ظهرت إحداثيات مثل 5200000, 2840000 فهي على الأغلب EPSG:3857، ولا تضعها في حقل خط العرض والطول.
• توقع اسم عربي لكل مكان: الخريطة تعرض name:ar إذا وجد، وإلا الاسم الأصلي، فالحل في بيانات OpenStreetMap وليس في الكود.

تدريب عملي

• احسب عنوان المربع لمدينتك في المستويات 5 و10 و15 بالملف tiles.py، ثم تأكد من النتيجة بالصفحة tiles.html بعد أن تقص مستطيل مدينتك.
• أضف إلى geo.py مسافة بين مدينتين في شمال أوروبا، مثل أوسلو وستوكهولم، وقارن خطأ الحساب على Web Mercator بما رأيناه في المملكة.

الخلاصة

وصلنا لنهاية الموضوع، وأهم ما فيه:

• خط العرض وخط الطول زاويتان على نموذج WGS84، أي EPSG:4326، وأغلب المكتبات تكتب الطول أولاً.
• Web Mercator، أي EPSG:3857، يجعل العالم مربعاً سهل التقسيم ويحفظ الأشكال المحلية، ويكبر المساحات بمعامل 1/cos²(lat)، أي 1.21 في الرياض و10.47 عند جرينلاند.
• في المستوى z يقسم العالم إلى 2^z × 2^z مربعاً، والمعادلة تعطي x وy، وجزؤهما الصحيح عنوان المربع وكسرهما مكان النقطة داخله، والرياض في المربع 12/2579/1758.
• المربع المتجه يحمل الطبقات: الأرض والماء والطرق والمباني والأماكن ونقاط الاهتمام، وكل اسم فيه بعدة لغات، ومربع وسط الرياض 99 كيلوبايت فقط.
• كل رقم بعد الفاصلة يقسم المسافة على عشرة، وستة أرقام تكفي أغلب التطبيقات.
• Haversine أخطأ بأقل من 0.1% بين الرياض وجدة، والحساب المسطح بأكثر من 7%، وعلى Web Mercator بقرابة 9%.
• GeoHash وH3 يحولان النقطة إلى خلية لها اسم، للبحث عن القريب وتجميع البيانات حسب المنطقة.

والخطوة التالية أن تعرض هذه المربعات من سيرفرك، ثم أن تحول العناوين إلى نقاط وتحسب المسارات بينها، وهذا ما في المقالات التالية من المرحلة.

دليلخرائطك على سيرفرك مع OpenStreetMap وMapLibreكيف تعرض خريطة تفاعلية لمدينتك من سيرفرك ببيانات OpenStreetMap، بالعربية والإنجليزية ودون مفتاح API أو مزود خارجي.دليلالبحث عن العناوين Geocoding مع Nominatimكيف تحول اسم المكان إلى إحداثيات والإحداثيات إلى عنوان على سيرفرك مع Nominatim، وما الذي يجيده بالعربية وما الذي يفوته.دليلحساب أقصر طريق وزمن الوصول مع OSRM وValhallaكيف يحسب محرك المسارات أقصر طريق وزمن الوصول على شوارع مدينتك، وكيف تشغله على سيرفرك وتقيس سرعته بنفسك.دليلكيف تبني تطبيق خرائط وملاحة مثل هدهد وGoogle Maps؟ما الذي يقف خلف كل ميزة في تطبيقات الخرائط، وكيف تبني نسخة أولى منها بأدوات مفتوحة المصدر وبيانات مرخصة بشكل صحيح.دليلالأنظمة العددية للمبرمجين: من Binary إلى Floating Pointلماذا يظهر القرص 465 بدل 500، ويصبح العداد سالباً فجأة، ولا يساوي 0.1 + 0.2 القيمة 0.3؟ الإجابة في طريقة تخزين الأرقام.دليلPostgreSQL لكل شيء: متى يكفيك ومتى تحتاج غيره؟كم عملاً يستطيع PostgreSQL واحد أن يقوم به في شركتك الناشئة قبل أن تحتاج Redis أو Elasticsearch، وما الإشارة التي تقول لك إن الوقت قد حان.

سجل التحديثات

• أكتوبر 2026: كتابة المقال، وتنفيذ التجربة بلغة Python 3.13 مع geographiclib 2.0 وmapbox-vector-tile 2.2.0 وh3 4.3.1 وgeohash2 1.1، وgo-pmtiles 1.31.2 على بناء Protomaps بتاريخ 20261008، وMapLibre GL JS 6.13.0 وnginx 1.29.8، وDocker Engine 29.8.0 وDocker Compose v5.5.1 على Windows 11 مع Docker Desktop.